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“The fundamental laws necessary for the mathematical treatment of a 
large part of physics and the whole of chemistry are thus completely known, 
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Chapter 1 

Introduction 


During the last decades, the topic of cluster physics became a genuine scien¬ 
tific and interdisciplinary field of research. This evolution was driven on one 
side by the fact that clusters provide a tool for tuning phenomena at nano 
or even atomic scale. On the other hand there was an academic boost in re¬ 
search due to evolution of experimental techniques which allowed for more 
precise experiments at nano-scale level. But perhaps the most important 
boost was given by the computer development. The fact is that since the 
invention of the transistor, the Moore’s law [1] has never stopped working 
and, in the 80's, the computers were already capable of dealing with quite 
difficult and demanding numerical problems. That was the period when the 
exponential era of many-body quantum theories began with applications like 
Density Functional Theory (which was basically dormant since 1965 when 
it was designed as practical tool) and mankind was finally able to simulate 
electron dynamics in the frame of different approximations in solid and in 
molecular physics. 

Having the numerical tool at hand, the theoreticians became interested 
in these wonderful objects, the atomic clusters, and in the 90's, one could 
already find in literature a considerable number of studies and reviews both 
on experimental and theoretical methods in cluster physics. 

The metal clusters took immediately the lead due to their simple struc¬ 
ture and the interesting features exhibited in the optical spectra. New be¬ 
haviour, not present in solid state physics, nor in molecular physics was 
seen and the interest in the interaction with electromagnetic fields ascended 
quickly in the hierarchy of subjects. 

Parallel with clusters another field of interest in physics gain considerable 
momentum in the same period: the lasers. This was driven by the Chirped 
Pulse Amplification (CPA) technique which was invented by Gerard Mourou 
and Donna Strickland at the University of Rochester in the mid 80's [2]. 
CPA putted the lasers on a track in which their intensity (until than almost 
stagnating) become exponentially growing in time and now, we are looking 
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forward to a generation of Peta-Watt lasers [3] . 

Naturally, lasers entered almost instantaneously in cluster physics as an 
easily tunable experimental tool for optical studies from which a large variety 
of spectra could be obtained and interpreted in terms of cluster properties. 
The theoreticians took their job seriously and soon, large codes implement¬ 
ing various theories were putted at work to reproduce the experimental 
results. 

As cluster and laser physics developed together into this beautiful flower 
of science, we are now at an edge where there is so many to know about 
what it has been done, but there are also serious theoretical questions about 
what and how to investigate further. 

It is the author’s opinion that a new view on the general many-body 
problem and in particular on the quantum many-body systems is needed. 
This could translate in an active research topic in conceptual and numeri¬ 
cal aspects but also in reaching new frontiers for laser-cluster phenomena, 
perhaps with many applications. 

While this might be true, the purpose of the present work is not to give 
such a new insight, but rather to establish a good, structured picture of the 
existing work with its pure theoretical or numerical aspects and to connect 
them with the experimental data. 

In literature, there is a series of reviews and books that discuss different 
aspects of the topic, but comprehensive studies on the theoretical approaches 
are rare. Usually, one has to come with a strong theoretical background or 
to span many different works to have a clear picture about the theory that is 
to be known. The truth is that the subject is so complex that not they, nor 
I can construct a complete presentation in such narrow spaces as reviews or 
academic theses. Perhaps, a series of books would cover it, but by the time 
one would write a volume, another volume should be written to cover what 
has been done in the mean time. 

Therefore, in the present work, I will try to fix a logical hierarchy of 
the theoretical methods to be used, on how are they connected with various 
experimental quantities or with different laser regimes and to provide enough 
further references that could allow a logical follow up into the subject. 

1.1 Problem setting 

In a standard thesis, the author can usually write a few pages on the central 
problem to be discussed. In my case, it is a challenge to point out what part 
of the theoretical methods used in laser-cluster interaction phenomena are 
more important and should be the core of the present work. It is rather that 
the whole concept of theoretical work in such complex systems is important 
and this is what a virtual reader should focus on. 

Without diminishing in any way the importance of the experimental- 


ist, in my opinion, the theoretician has a special challenge to overcome in 
the sense that is mandatory to have a wide spectrum of ” interdisciplinary” 
knowledge. For our topic, a moderate mathematical background is needed 
in terms of analysis, algebra, geometry, etc. This should be accompanied 
with a good understanding of the electro-magnetism, at least in the classical 
sense, since almost all interactions are of electro-magnetic nature. On the 
other hand, even if a lot of theoretical methods employ, for the sake of sim¬ 
plicity, classical or semi-classical methods, the presence of quantum effect 
is pregnant in cluster physics, and so, the need for a strong background in 
quantum statistics, is stringent. This quantum methods, superpose with the 
quantum chemistry and the theories from mesoscopic solid state physics in 
a strange way, but still, the topic remains self consistent. Not the least, a 
cluster can be viewed as a quantum plasma, especially at high temperatures, 
as it is the case in strong lasers fields, therefore, a good knowledge of plasma 
physics and in particular quantum plasmas is required, given the fact that 
many of the theoretical methods overlap, at least conceptually. 

As it is hard to say which theoretical method is more important, it is also 
hard to extract from the field of laser-cluster interaction a particular matter, 
or a goal to be achieved above all others. Therefore, I stress again that the 
main purpose of this thesis is not to give insight (as usually is done) in some 
narrow problem presenting the achievements of the author. In contrast, it 
was designed to be a pastel of how the zoo of theories grows on the physical 
problem and reflects how little the author has achieved while studying this 
matter. 

Finally, if there would be to mention a certain point in which the existing 
theory slips, is the fact that in strong laser fields, the scales involved span 
a few orders of magnitude, fact which makes any attempt to study the 
system with quantum methods an impossible task, the last resort measure 
being the classical methods which miss essential quantum phenomena. The 
compromise has to be done between classical and quantum theories and a 
mid-path is, in my opinion, the main problem in the present state of the 
research. 

Writing this thesis required a lot of reading, but beside all the references 
cited during the following pages (and many others) it is important to men¬ 
tion some fundamental books and articles that stood at my bedside during 
the last two years m El El 0 Eli and some related (well written) PhD thesis 

uni EH EH [i3i m ES] 

1.2 Applications 

For the author, the most important aspect of the field of laser-cluster cou¬ 
pling is that it represents on of the richest systems in which the problem 
of radiation-matter interaction can be studied from a quantum many body 


perspective. 

But the truth is that the scientific interest in the subject was just par¬ 
tially driven by this academic interest. The applicative possibilities were 
the ones that drove strongly the field. Just for a visual impact we should 
note some references which discuss the applications of clusters in medicine 
dSIITlIIHlIIlEolEI], chemistry [22l [23l [M] , optics [25l [26] , microelectronic 
1221, etc. The list can go on, but it is less focused on the laser part and more 
on how sole clusters could be used. To have a more detailed enumeration of 
what can be achieved (until now) from clusters irradiated with strong lasers, 
I should mention: 

1. Production of shaped ion energy spectra [28l (29] 

2. Production of highly energetic charged ions jJO] . 

3. Ejection of hot electrons m 

4. Emission of extreme UV and X-ray photons [aiiiMllMj 

5. Ultrafast X-ray diagnostics [35l [36] 

6. application of cluster explosion in inertial fusion confinement 1371 . 

mm 

7. Avalanches of ionized electrons can be used as a mechanism of damage 
in solids with applications in material processing, microfluidic devices, 
nano-surgery [lO] 

8. A novel application which until now has been used only in small 
molecules m would be to use the Coulomb explosion imaging ap¬ 
proach to nano scale systems. 

9. High harmonic generation |42[ I43j . 

More detailed discussion and general situations can be found in |44| . 

1.3 Thesis structure 

Beside the present Introduction [T| this thesis is structured in 3 chapters. 

The second chapter is designed to give a short definition with appro¬ 
priate pictures to the concept of the cluster. The picture is complemented 
by a shorter description on how laser pulses are modelled in the domain of 
interest. A gross classification of the main dynamical regimes of interac¬ 
tion is given with its quantities to be connected with the experimental and 
theoretical methods later on. 


The third chapter is the heart of this work. It gives us the structured 
picture of the existing theories used in cluster physics starting from first prin¬ 
ciple and axiomatic Quantum Statistics (Liouville-Von-Neumann equation) 
and going through different levels of approximation until the semi-classical 
and even classical methods are obtained. For logical reasons, a short, formal 
derivation in cascade is tried to be exposed in there. Since this is just a mas¬ 
ter thesis, one cannot hope for refined details about any of the presented 
theories (about many of them books have been written) but rather for a 
schematic image which should allow one to choose easily the appropriate 
method to be used. 

Being a theoretician is hard... Your experiment is the numerical simula¬ 
tion. For this reason the fourth chapter [^exists. It is a representation of the 
numerical work (which, as always, took almost 80% of the time devoted to 
the thesis) that has been done in the process of understanding and practising 
the theory. It was the author’s intention to apply for different phenomena or 
regimes of study in cluster-laser physics as many theories as possible and to 
give insight about the numerical methods involved. If someone would read 
this work, that is the place where fancy coloured pictures are! 

The last chapter 4.4 provides a set of conclusions, future plans and per¬ 
spective of the problems discussed in the thesis. 



Chapter 2 


Clusters & lasers 


2.1 Clusters: the world of not too few, nor many 
enough 

Before entering in the technical details of laser-cluster interaction and the 
associated theoretical methods, it is necessary to have some understanding 
about what is a cluster in general and what are the characteristic properties 
to be expected in such systems. 

Clusters are, by definition, aggregates of atoms or molecules with regular 
and arbitrarily scalable repetition of a basic building block. Their size is 
intermediate between atoms and bulk. One could thus loosely characterize 
them with a formula of the form where (3 < n < 10®“'^). [3] 

In a pedestrian view, one could imagine that a cluster is some collection 
of atoms (many times the same type of atom) that stick together in some 
random shape. First question would be: why not to call them molecules? In 
the end, a molecule is also a (small) collection of atoms. Well, the difference 
is done by the larger number of atoms in a cluster and the properties therein. 
Perhaps the most pregnant difference is that the usual molecules has a small 
number of isomers and a precise (easy to identify by numerical means) con¬ 
figuration for the ground state. In cluster, this is no longer the case since 
even small clusters can present large number of isomers. For example Aris 
it is known to have hundreds of isomers. On the other hand, although the 
quantum effects are present (and essential) in cluster dynamics, many times 
they behave in a coherent manner giving rise to collective phenomena, less 
present in molecules. 

At the other extreme, one could ask: why not to consider very large clus¬ 
ters as being solid bulk and study them with the associate methods? Well, 
again it is a matter of quantum effects. It can happen that, even large, the 
arising quantum phenomena to cannot be neglected in a classical manner. 
But the most important difference between usual solid state systems (as 
crystalline lattices) and clusters (which are finite systems), is that the latter 
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Figure 2.1: A pedestrian view on a gold cluster 


lacks the periodicity. There might be symmetries as in spherical clusters or 
fullerenes, but the translational symmetry (the usual solid state periodicity) 
is generally not present, therefore, it is hard to speak about band structure, 
etc. 

Regarding the cluster material, as said before, it is standard to study 
clusters of a single type of atom. Historically, the first to emerge and by far 
the most studied are the metal clusters {Na, Li, etc.). This happened for 
many reasons: the quasi-free electrons, the easiness of producing them, the 
strong optical response, etc. Another class of clusters that attracted great 
interest during time is the class of fullerenes which are carbon compounds 
with high symmetries and astonishing physical and chemical properties. In 
the matter of intense laser-cluster interaction, has become customary in the 
past years to study the response of rare gas elements clusters {Ar, Xe, Ne, 
etc. ) especially due to the fact that they are a source of many exotic 
phenomena (the X-ray production, for example). 

The theoretical approaches which will be discussed in Chapter will be 
used further in the numerical chapter]^ of this thesis to investigate properties 
of some types of clusters, in particular the Nan, Cn and Xe„. For this 
reason, I shall focus next on describing the main physical scales for this 
three cases. 

2.1.1 Characteristic scales 

In general is important to have a good intuitive sense about the magnitude 
of the observables in a physical system since it can give you insight about 
the applicability of various theories. In this matter, we should describe in 
this subsection some specific quantities characteristic to clusters. Similar 
discussions can be found in [Zlilll] 


1. The number of constituents. The number of atoms in a cluster 
spans almost 6 orders of magnitudes. With this, problems or simplifi¬ 
cation arise. Generally, as in nuclear physics, there are specific clusters 
(the so called magic clusters) which are more abundant in universe and 
this can be seen also from mass spectrometry experiments. The exis¬ 
tence of magic clusters is a consequence of the quantum nature of the 
electronic structure which allows for the presence of quantum shells. 

2. Wigner-Seitz radius r^. The concept is used as the radius of a 
virtual sphere in which an atom from some material is contained in 
the bulk regime. Therefore, it is related with the type of material and 
has a weak physical meaning for (small) clusters. Still, many quantities 
can be expressed or normalized in terms of it. 

3. The size. As one can imagine, when you are dealing with a few atom 

cluster, the specific dimensions are lA. But as the number of atoms 
increases, one can hnd at the superior limit, clusters with a diameter 
of ~ lO^nm. For spherical shaped clusters, the radius can be roughly 
expressed as Rq = where N is the number of atoms. From 

there a gross electronic density can be computed as po ~ 3/(47rrg). 

4. The time scale. The time scale is very important in the dynamical 
regime. In metal clusters for example, where the electrons are quasi 
free and they respond very easy (and fast) to external fields, the dy¬ 
namics has a specific time of ~ 1/s. This approximate value is a gen¬ 
eral magnitude order for electrons in all kinds of clusters. This aspect 
reflects in experiments and numerical simulations, where the electron 
dynamics is usually investigated for at most a few hundreds of /s. On 
the other hand, the numerical time steps needed are around 10“^/s. 
The implications of this values are important for the laser-cluster in¬ 
teractions because it basically sets a specific range of frequencies at 
which the optical response of the cluster is resonating. As we will see 
in Sect. |2.3t the dynamics of ions can be neglected for weak lasers 
(smaller then I < 10^^VF/cm^). But when the power exceeds such 
values, the response of ions to the external fields is no longer negligi¬ 
ble and we must deal with different ionic motions, from vibrations to 
escaping particles (as in a Coulomb explosion). Nonetheless, having 
in mind the ratio between the mass of a nucleon and an electron, it is 
clear that the ions have larger characteristic time scales with a generic 
value of tens of fs or even ps. 

5. Resonance frequency. A significant part of the theoretical and ex¬ 
perimental investigations of clusters is concerned with the optical or 
more generally, the dynamical response (linear or non-linear). The 
purpose is to find maxima in the response spectrum which mark the 


energetic resonances of the system. The general case requires compli¬ 
cated approaches, but a first glance on the magnitude of the resonance 
frequency (related to Mie plasmon in metal cluster) can be obtained 
modelling the cluster as a spherical metal drop and by Mie’s theory 

m- 
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6. Landau fragmentation/collisions. In clusters, this phenomenon, 
known from classical plasma physics is related with the coupling be¬ 
tween plasmons (collective) and single-particle excitations with en¬ 
ergies close to hojMie- Roughly, the characteristic time of Landau 
damping can be expressed [l6l 07] tl = /vp, where vp = 

(h/m)(97r/4)^/^ is the normalized Fermi velocity. 


7. Electron-ion collisions. Is a phenomenon which appears due to non 
zero temperatures in clusters and scales as Tei oc T~^ at low temper¬ 
ature, due to electron-phonon scattering [JH], while at high tempera¬ 
tures r oc [19| . 

8. Electron-electron collisions. Is a phenomenon related with the 
termalization of a cluster, or the evaporation of electrons. At law 
temperatures the characteristic time goes like Tee cc T~^ accordingly 
with Fermi liquid theory [50|. 

9. Ionization potentials. This quantity is in general related with the 
energy of HOMO (the highest occupied molecular orbital) and it has 
lower values in metals. Still in general clusters, a wide range of val¬ 
ues can be found, an roughly approximation being the order of eV. 
An interesting aspect is that the IP in a cluster can be dramatically 
different from the same IP in a single atom from the same element. 


10. Critical laser intensity. It represents the needed intensity of the 
laser field to induce ionization in a cluster. Even though tunnelling 
ionizations can appear, the quantity is closely related to the IP. 

11. Thermal de Broigle wavelength \p. Stays as a measure of the 
quantum effects in a thermal plasma. In an intuitive representation it 
can be seen as the spatial extension of a particle. When this dimen¬ 
sion of the plasma increases and Xp becomes less than interparticle 
distance, then the quantum effects fade away. 
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Xb = 

mvT 

= ^ ( 2 . 2 ) 

V m 

12. Fermi temperature. The quantum statistics of free particles (Fermi 
Dirac) allows for a direct link between the density and the fermi energy, 
therefore, one can define a quasi-temperature similar with the classical 
case of ideal gas; 
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13. Regime parameters. We can define the following quantum coupling 
parameters [1]: 
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These two non-dimensional parameters define our regimes and con¬ 
sequently the theoretical methods to be used: y <C 1 —)• classical, 
X 3> 1 —)• quantum, qq 1 ^ collisional, x 3> 1 —collisionless. 


As a short picture about the values involved, in Tabel 2.1.1 
some values of some quantities vs. specific elements. 


one can see 


Element 

IP[eE] 

rs[M 

^Mie[ey] 

Icrit[W/ cm?] 

Na 

5.14 

2.1 

2.8 

3 X 10^2 

Ag 

7.58 

1.59 

0 

1 X 10^3 

C 

11.26 

1.21 

20 

6 X 10^3 


2.2 Lasers 

From the perspective of laser-cluster interaction, we are less interested in 
how the laser field is produced, but rather in what ’’comes out” from it. 
Now, everybody knows that the laser field is an electromagnetic (EM) wave 
in its most pure definition but there are various parameters that can, or 
cannot be, neglected when speaking about its effects on molecular scale. 

















First of all, the laser field is completely described by the knowledge of 
the two vectorial fields E electric and B magnetic which obey the Maxwell’s 
equation. Second, a rough length scale of Inm is much less than the spatial 
extent of the laser profile (which is usually gaussian in section), therefore, the 
dipolar approximation holds: the laser field is considered constant over the 
region of interest. On the other hand, the overall angular momentum and 
charge current in a cluster is small and consequently, its coupling with the 
magnetic component of the EM wave can be neglected. This considerations 
lead to a natural modelling of the laser in only one constant electric field: 

£{r, t) = £oezf{t)cos{u;iast + (2.6) 

Now, it is clear that is the versor of the Oz direction in space which is 
chosen arbitrary just for a simplification of calculus (even though it maps to 
the linear polarization), £q is the magnitude of the electric field (a constant) 
while the rest of the formula |2.6| represents just time dependent terms. The 
cos term models the oscillatory behaviour of the EM wave which has as a 
first parameter the laser frequency ojias which relates to the energy of the 
photon carrier by the trivial hwias- The term is a phase term which can 
be written roughly as (/)(t) = (f)ce + (/5/2)t^ + {'y/3)t^ + C>(t^) with (/)ce the 
carrier-envelope phase, j3 ,7 linear and quadratic chirps. 

The instantaneous intensity of the laser field is easily expressed I{t) = 
ceoTo/2/(t)^. The function f{t) models the temporal shape of the pulse. A 
common form for it is the Gaussian field profile f{t) = exp(—(t/r)^). 

Now, as we shall see in Chapter many quantum approaches on the 
electron dynamics do not work with the electric force acting on them, but 
with the effective potential in which the electron’s wave function behaves. 
Therefore, it is useful to have an associated electric potential for the laser, 
and in the frame of dipolar approximation it can be expressed as: 


Vias{r,t) = e£r 


As a key parameter for the classification of coupling regimes between 
lasers and electrons is the ponderomotive potential which reads: 


Up = 




= 9.33 * 10"^^/o[lT/cm^](A[^m])^eE 


This quantity is a measure of the relativistic effects to be expected. 
When Up becomes comparable with the rest energy of the electron, then the 
relativistic regime becomes active. In general, through this thesis, I do not 
refer to such extreme cases. 

It is worth mentioning that a more general description of the laser field 


is the modelling of a train of pulses [5T] and the relation 2.6 generalize itself 
to : 




£{f, t) = f’oe^ ^ g{t)cos^{^ —^7r)0((t - t«)(T - t + ta))sin{ujiast\2.7) 

CX.=0 

ta = Ar + aTtrairi2.8) 

g{t) = exp[-4ln2- —](2.9) 

r = ATt^i^X^.lO) 


2.3 Laser-Cluster coupling regimes 

It is important to have some understanding about what different laser pa¬ 
rameters bring on the table in terms of phenomena in clusters in order 
to have later a good sense about the theoretical method to be applied. 
The key feature that spans all the regimes to be discussed below is the 
(photo)ionization. This subsection was highly influenced by [7]. 

From an atomic perspective, the photo-ionization can be divided in two 
categories: vertical ionization and optical field ionization (OFI). The vertical 
ionization consists of a basic single or multiple absorption of photons by the 
electrons from the atom, and consequently, their transition on higher bound 
or free states. In general, a process of multiple photon ionization (MPI) 
which involves n photons has a specific reaction rate that can be expressed 
as F = where Cu is the absorption cross-section and I is the intensity of 
the laser field. On the other side, external fields with a slow time oscillation 
can be treated as static and so, the deformation of the atomic potential has 
a long enough life to allow ionization through tunnelling of the energetic 
barrier. This is known as OFI. An important (hystorical) parameter here 
is the Keldysh parameter [52] 7 = Eip/2Up which indicates the type of 
process that dominates: 7 ;§> 1 means MPI and 7 < 1 means OFI. 

When moving to clusters, the fact that in the system are more atoms 
makes the simple picture of atomic potential more involved. The interaction 
between atomic electrons and their overlapping deforms the effective poten¬ 
tial in which they move and the details of ionic and electronic structure 
become crucial. The simple fact that two atoms are bind together, usually 
means that the potential is lowered between them and the effect known as 
charge-resonance-enhanced ionization (CREI) arises. This comes in package 
with an increase in the ionization rate. 

Essentially, in clusters, there are two types of ionization: inner ioniza¬ 
tion and outer ionization. Inner ionization consists of vertical absorption of 
photons which excite electrons from bound states (core electrons) to valence 
states but still bounded to the cluster (quasi-free electrons). The outer ion¬ 
ization goes a step further and represents the transition of valence electrons 
into the continuum (free spectrum). The main channel of energy absorp- 






Figure 2.2: A schematic representation of the effective microscopic potential 
yRS (blue), the deformation induced by large external fields Viaser (gray, 
dashed) and the main types of states: core electrons, valence electrons and 
continuum states; In green, occupied states/ in black, holes 


tion is the direct transfer from laser field to the electrons. Still, during the 
dynamics of the electron cloud, the self-consistent Coulomb potential may 
play itself a role in the single particle excitation, creating metastable states 
in the (quasi) free spectrum. 

Following the laser parameters (in principle the intensity) we can draw 
three different regimes of interaction. 

2.3.1 Linear regime 

Is characterized by laser fields with small intensities, usually bellow 10®IF/cm^ 
Correspondingly, the Keldysh parameter is very large since the pondero- 
motive potential is much lower than the ionization potential 7 ^ 1 . The 
processes are mainly frequency dependent. From them we distinct the photo¬ 
absorption related with the quantity called optical response (cross section) 
and the single photo-ionization related to the photo-electron spectroscopy 
(PES). From a theoretical point of view, this regime can be tackled with 
linearized methods. 

The optical spectra is constructed from the dielectric function e{oj) which 
in turn is proportional with the Fourier Transform of the total dipolar mo¬ 
ment. A more suitable quantity to represent for non-linear features is the 
power spectrum Siuj) [53] : 


















( 2 . 11 ) 

( 2 . 12 ) 


cr{u!) oc a;3'(e(w)) oc ^^[^^(w)] 

S{uj) oc |Zi(a;)|^ 

As we shall see in Chapter the PES can be computed from single 
particle methods recording the time evolution of the orbitals at a ’’far-away 
point”. From a physical point of view, a PES spectrum is a picture of the 
density of states in the cluster. The main problem with the PES is that 
experimentally is hard to achieve for core electrons. But if one is interested 
in the single photon processes and linear regime, these types of excitations 
are not possible. 

More insight in the geometric geometrical structure and the active or¬ 
bitals during the dynamics can be achieved in photo-emission spectra, from 
a spatial analysis. More precisely, it is investigated what is called Photo¬ 
electron Angular Distribution (PAD). This can be quantified in the differ¬ 
ential cross-section which is expanded in a series of Legendre polynomials: 

^ = ^(1 + P2P2{cos{ 0)) + p^Pi{cos{e) + ...)) (2.13) 

ail 47 r 

At the first glance, it can be seen that the /32 parameter is a reflection 
of the orientation of the emission over the polarization of the laser field. 
A /32 = 2 means emission parallel with the laser electric field, /32 = 0 an 
isotropic emission and /32 = — 1 emission perpendicular on the polarization. 
An important issue with the PAD is that is photon dependent (in the sense 
of frequency) [M] . 

While PES and PAD give us information about the ground state struc¬ 
tural properties of the cluster, a final tool in the study of single photon-linear 
regimes can be developed to investigate the dynamical features, namely the 
Time Resolved Photoelectron Spectra (TRPES). 

2.3.2 Multiphoton regime 

Is characterized by intensities ranging in the 10® < I < 10^3IT/cm^. In 
terms of Keldysh parameter 7 < 1. This time, all the characteristic pa¬ 
rameters of the laser field (w, I, E, f{t)) become equally important in the 
dynamics. The MPI takes the leading role in the ionization and as a conse¬ 
quence, the phenomena of second harmonic generation appears. The optical 
spectrum presents non-linear features. All the linearized theoretical ap¬ 
proaches break at these intensities and full propagation schemes (Vlasov, 
TDDFT, TDHF, TDTF, etc.) must be employed. 

When in this regime, the multiphoton processes are always accompanied 
by single-photon ones. In principle from PES experiments (or calculations) 
one can extract the single particle energy by the relation Ekin = i^huj -|- eo- 
Reversely, if one knows eo and the laser frequency, the type of phenomena 


single/multi can be recovered computing v. By default, laser fields with fre¬ 
quencies bellow the ionization threshold can ionize the cluster only through 
OFI which has a very small probability. 

A very robust phenomena captured by photo-ionization spectra is the 
plasmon. More details about it will be given in|^ but in principle, this is a 
collective phenomena characteristic to metal clusters. It remains a pregnant 
resonance even in non-linear regimes where the single particle features can be 
washed out by the laser intensity. Being a collective phenomena it represents 
a gate for the resonance and a good energy absorption of energy in the 
cluster. Various types of damping behaviours can fragment the plasmon 
and give a spectral picture with many peaks around uJmie- 

2.3.3 Strong field domain 

Is achieved in the 10^^ < I < jerr? range of intensity. Although there 

are studies which invoked the motion of ions during the dynamics, in general 
it is considered that the ions have a small amplitude and slow motion for 
I < lO^^VF/cm^. When the laser exceeds this limit, the ionic background 
cannot be considered anymore frozen and its dynamics must be taken into 
account through theoretical methods related to the Molecular Dynamics. 

The guilt for these violent dissociation of the cluster is a very efficient 
energy absorption. For intensities around lO^^VF/cm^ the mean energy per 
atom can be of ~ 100/ceV. This energy absorption ionizes very rapidly the 
cluster stripping the valence electrons and giving some boost to the inner 
electrons and ions. Having a considerable positive charge, the cluster dezin- 
tegrates itself even if the laser pulse was short due to the Coulomb repulsion. 
This self-consistent interaction gives in turn the so called Coulomb explo¬ 
sion. In turn, the hot quantum plasma created ejects with high velocities 
electrons, ions and photons (from recombinations). A schematic picture of 
the cluster explosion is presented in |2.3[ 

From the ejected ions, a first feature is the presence of highly charged 
ions. Different studies have reported even ions with a ~ -|-20 charge. Be¬ 
side high charge, the ions come in package with high energies at the level 
at even IMeV [301. As amazing as it is, these energies have driven imedi- 
atelly the researchers to study the possibility of using them in nuclear fusion 
applications. 

Regarding the electron dynamics, their inner ionization into a nano¬ 
plasma gives the possibility of recombination. In turn, the recombination of 
electrons with ions with empty core shells leads to the emission of photons, 
in particular, energetic photons. So, X-ray production becomes an active 
possibility together with high-order-harmonic generation. 



Figure 2.3: The cluster explosion is represented in a schematic view. The 
incoming laser pulse (in red) excites the ground state cluster (with yellow- 
electrons, blue-ions) which expels electrons and ions during the dynamics 












Chapter 3 

Theoretical methods 


The topic of laser-cluster interaction requires a wide set of knowledge from 
various fields of Physics. Nonetheless, the main requirement is the knowledge 
of Quantum Physics, coupled with some good understanding of the classical 
Electrodynamics. The present chapter, being a logical (hierarchical) picture 
of the existing theoretical methods used in cluster laser-physics, will follow a 
line which start with the fundamentals of Quantum Physics, more precisely, 
the basics of Quantum Statistics. In the end, the so-called classical methods 
will be discussed, but as the final point of a line of approximations. 

Almost every many-body theory that will be discussed can be derived 
in more than a single way, following different representations of quantum 
mechanics, or different mathematical approaches. Still, the practical end 
is always the same, so we will stick with a self consistent presentation, in 
parallel, a sufficient number of further references being given. 

Regarding the other realm of modern physics, namely the theory of rel¬ 
ativity, we will not refer in general to this kind of effects since they be¬ 
come important only in the extreme ultra-intense lasers, with intensities 
above 10^®IT/cm^. Still, almost any model which will be discussed can 
support relativistic versions or corrections without modifying the main idea 
of the theory. For example, methods as Density Functional Theory which 
are treated with Schroedinger-like equations will, naturally, be extended to 
Dirac-like equations. Furthermore, the QED effects will be also neglected 
during the presentation. This not just a choice of the author, but actively, 
this kind of studies in clusters are almost singular since the treatment is too 
complex, while the effects are negligible Also, while the main part of the 
dynamics in a cluster is taken upon the electrons, the spin will be taken into 
account only where is absolutely necessary. Otherwise, the extension of the 
theory to spin is almost trivial. 

For all these reasons, this chapter will start with a discussion about the 
fundamentals of Quantum Statistical Mechanics (QS). 
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3.1 Quantum Liouville-von-Neumann equation 


Any science is defined by its set of postulates and its mathematical appara¬ 
tus, therefore I shall start with the axiomatic logic of QS before treating any 
specific theory. Moreover, as we shall see later, even in the strongest exter¬ 
nal fields, a separation of the cluster will be at hand in nuclei (possibly ions) 
which are subject to (quasi) classical behaviour and electrons, subject to a 
more or less sophisticated quantal treatment. As the complex treatment is 
required by electrons, it is natural to discuss only theoretical methods that 
deal with systems of identical particles, in particular, fermions. 

From a historical perspective, the hrst rigorous axiomatic of quantum 
mechanics was established by P.M. Dirac |55L [56] and it is a viable con¬ 
struction for the case of fully known, pure states. Unfortunately, during 
dynamics (especially strong dynamics) we deal with finite (possible large) 
temperatures, non-equilibrium phenomena in which the states of the system 
are no longer pure. Therefore, a more natural axiomatic, which should in¬ 
corporate statistical aspects of the dynamics, is needed. This was brought 
on solid mathematical grounds by J. von |57!. 

While the introductory quantum mechanics associated with pure states 
of particles works with the concept of wave-function or the more abstract 
notion of vector in an Hilbert space, the physics of mixed states is described 
in terms of density operators, which will be denoted through the thesis with 
T. There is no need or purpose to discuss the explicit axioms in here but 
some words should be said about the properties of the density operators and 
their dynamics. 

First of all, F is defined on the Hilbert space % of the system (as are 
the vectors for the pure states) and there are some characteristics that any 
density operator must obey. First, the trace must be 1, Tr(r) = 1, or more 
generally, N where N is the number of particles from our system. Second, 
in order to have real observables associated with our system, F must be 
self-adjoint F = F^. Later, we will see how there are interpretations of F 
in terms of density of probability and for that reasons the operator F must 
be positively defined and bounded. In principle, if one knows F , then, any 
observable A can be computed in terms of the associated operator A taking 
the trace in the Hilbert space of the system (A) = Tr{AT). 

Beyond these properties of F and the first axioms of QS which are ba¬ 
sically the same with those from the usual QM, it is important to start 
from the last of QS’s axioms, the one that describes the quantum dynamics, 
namely, the Liouville-von-Neumann equation: 

ihdtt=[H,t] (3.1) 

This is the quantum version of the classical Liouville equation. Both, 
are consequences of the Liouville theorem that states that in a Hamilto- 


nian system the distribution function is constant along the trajectories, or, 
equivalently that in the phase space the volume element is conserved. While 
the parallel between classical and quantum Liouville eq. is stringent from a 
visual level, it remains a single (apparent) difference in the fact that the clas¬ 
sical Poisson bracket is replaced with a commutator. This change is known 
as the correspondence principle and is a consequence of the quantization of 
the classical phase space in a Hilbert space. 

Now, let us denote for future purposes our system as having N iden¬ 
tical (indistinguishable) particles at the coordinates xtv = {xi, X 2 , ■■■, x^} 
described by a field operator T(xAr). Having this, the density operator can 
be formally written: 


f(xAr,x)v) = T(xAr)Tt(x)v) (3.2) 

Another concept that soon will be of great use is the reduced density 
operator of order ”s” which is nothing else than a partial trace from the 
entire P, on the sub Hilbert space of N — s particles P^ = Tr^^^P. If one 
prefers a representation of P in the coordinate space, than, it can write 
rs(xs,x's) (which in literature is found as distribution kernel ofT) as: 

fs(xs,x'^) = ^ ^ ^ y T{xN,x^)dx'g_^_i...dx'^dxs+i...dxN (3.3) 

Now, let us assume that our hamiltonian operator describes a system 
with at most two body interactions = ^il‘^xnVYl,i v{xi)+l/2 
Xj\) (this is always the case for electrons in clusters). Taking this partial 
trace of the Liouville eq. |3.1| we obtain an infinite hierarchy of equations, 
the so called BBGKY (Bogoliubov-Born-Green-Kirkwood-Yvon) hierarchy 
analogous to the one obtained in classical statistical mechanics : 

= [.H^f,] + Q^(f,+l) (3.4) 

■ 5-2 ^ 

= + Xj) (3.5) 

i=k k^j 

s 

Q^ts+i) = {N-s)Y, Trs+i{[(f{xk - x.+i), f.+i]} (3.6) 

k 

More or less cumbersome, the system is exact but infinite and there is 
no way to solve it. On the other hand, there is no use of solving it, since 
the knowledge of the entire density matrix is obsolete in practice. We are 
usually interested in macroscopic quantities as density, current, multipolar 
moments, etc. and this kind of information is accessible just from the first 
order reduced density operator Pi. But one cannot solve the corresponding 


equation for fi due to the term which is intrinsically dependent on the 
all other higher order density matrices. 

The Liouville’s eq. 


3.1 allows for a variable number of particles during 


the dynamics. This aspect is preserved in the BBGKY hierarchy and it 
becomes obvious that one can have creation or particles, usually in strong 
external fields. Other important properties of BBGKY are that the system 
can be formally solved, due to linearity and it contains an intrinsic math¬ 
ematical time reversibility. Regarding the conservation of energy, unlike 
other kinetic approaches (Boltzmann, Landau, or Balescu [58] since they 
are derived under the assumption of zero three particle correlations), the 
system conserves the total energy. 

There is some mapping between the density operator and the single 


m 


particle representation of a statistical system through the relation 3.7 
which Ti is expressed as a superposition of projectors in the single particle 
Hilbert space, weighted by some probability coefficients {pj}- 


Ti = 


(3.7) 


The BBGKY system of equation is the basis for part of the methods 
which will be discussed further. In particular we will start from the equa¬ 
tion for r 1 and use some approximations to decouple the hierarchical depen¬ 
dence with higher order densities. Just to have a visual picture about what 
is to emerge from these approximations, in Fig. 3.1 there is represented an 


organizational chart with the main theories and their relation with higher 
theories. Starting from the Von Neumann axiomatic and the density op¬ 
erator formulation of QS we see how we have on one hand the quantum 
Liouville’s equation from which, passing trough the BBGKY hierarchy and 
some approximations two lower theories can be obtained: the Hartree-Fock 
(HF) theory and the Quantum Wigner (QW) equation. Both are at the 
same level, being equivalent conceptually but used under different represen¬ 
tations. On a parallel level is the Density Functional Theory (DFT) which, 
apparently, is not related with QLiouville. This fact is not true, but the path 
on which DFT is derived has no direct connection with the latter mentioned 
equation. Both HF and DFT, being in practice, single particle methods, al¬ 
low for a linearisation and a description of normal modes in terms of single 
particle wave functions and energies. From HF it is obtained the so called 
Random Phase Approximation (RPA) while from DFT the Linear Response 
DFT (LR-DFT), both having formally equivalent results. 

A different theory can be extracted from QW Eq., taking the semi clas¬ 
sical limit fi —7- 0 and retaining the zero or the first order in h, namely the 
Vlasov equation. This is equivalent with the classical Vlasov equation used 
in plasma physics, but quantum effects can be introduced through mean field 
potential (taken from DFT) or through collision operators that describe the 






Figure 3.1: An hierarchical map of the main theoretical methods used in 
cluster physics 


semi classical quantization of the phase space (Uheling-Ulenbeck is one type 
of such operator used to mimic the Pauli principle). 

Going even lower in the tree of approximations, one can derive a Quan¬ 
tum Hydrodynamic Model (QHM) from three different perspective: inte¬ 
grating the Vlasov equation on moments in the momentum space, using the 
abstract (functional) Euler-Lagrange equations of DFT or using a Madelung 
transform on the orbitals from the DFT/HF equations. The last level of ap¬ 
proximation is designed for the intense laser-cluster interaction and large 
clusters, namely the nano-plasma model or classical molecular dynamics 
(CMD). 
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Table 3.1: Schematic table with the gross applicability of the main theories 
to be discussed: Hartree-Fock 3.3, Density Functional Theory 3.4 Random 


Phase Approximation 3.7.2, Thomas Fermi (and Orbital Free extensions 


3.6), Vlasov 3.5.2[ Molecular Dynamics 3.2.2 and Rate Equations (from 
nano-plasma model) 3.8.2, E/N stands for excitation energies, N for the 


number of atoms in the cluster and I for the intensity of the laser field. 
L/M/S =Linear/Moderate/Strong regimes. 


The tabel 3.1 shows a schematic view of the applicability of each of 
these theories. As one can see, the large number of atoms coupled with 
large intensities has a poor representation. While there is no place or space 
to describe any of these theories in detail, all of them admit extensions and 
improvements which will be only mentioned with appropriate references. 


3.2 Quantum first principles in atomic systems 

Clusters, as atomic systems are quite complicated things to study. Think 
for a moment that you have an Xn cluster which will contain by default in 
a neutral state n nuclei and nZ electrons moving around accordingly with 
the quantum mechanical rules. 

As a first step in the investigation of dynamics (or ground state) we 
should write down a Hamiltonian H for our system. Accordingly with all 
we know from elementary classical and quantum mechanics the Hamilto¬ 
nian should contain a kinetic and a potential energy term, the later given 
by the Coulomb interaction between electronic and nuclear charges. Let 
us denote with {R/} = 94 the nuclear coordinates, with (r)} = r the elec¬ 
tron coordinates and with 'k(94, r) the total wave function associated with 
electron-nuclei system. With those, the Hamiltonian and the Schroedinger 
equation can be written as 












H^{D\,x) = E^{^R,x) 


(3.8) 

(3.9) 

(3.10) 


H — Te + Tjv + Vee + Vnn + Ven 
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(3.11) 


(3.12) 



ZiZj 


(3.13) 


(3.14) 


Where Mj is the mass of the I — th nuclei, m is the electron mass and Zj 
is the atomic number of the I — th nuclei. Now, our simple n atom cluster 
involves only for the ground state, to solve this Schroedinger [3^^ equation, 
which is a partial differential eigenvalue problem in 3n(l + Z) coordinates. 
Obviously, this is an impossible task both analytic or numeric. Therefore, 
stated as it is, the full atomic problem in a cluster is superfluous. The next 
section will present one of the basic approximation in any molecular, cluster 
or solid state physics problem, the Born-Oppenheimer approximation. 

3.2.1 Born-Oppenheimer approximation 

Now, having the |3.9| Hamiltonian written down and explained, we should 
see briefly what is the Born-Oppenheimer approximation |59] and how it 
can be derived. 

BO: In an atomic system, the electrons-nuclei problem can he separated 
in two distinct problems. Consequently, the factorization 'I’totai = tpeiectron x 
f^nuclei holds. 

Let us look first at the magnitudes of some constants that appear in the 
Hamiltonian. A nucleon is roughly 1800 times heavier than an electron, so 
we can safely state that Mj ^ m. This means that in a classical view of our 
system, the nuclei move much slower than the electrons (being subject to 
forces of the same magnitude. Coulomb forces, but heavier particles). For 
this reason we could, as a first step, neglect the kinetic energy of the nuclei 
and write the so called ’’clamped nuclei” Schroedinger equation: 









(3.15) 

(3.16) 


Hel — Te + Ven + Vee 
= Eel-lpeli^; 91 ) 

Where, this time, the 'I'ez(i^;91) wave function stands for the electrons 
and the nuclei coordinates are considered only as parameters appearing in 
the electron-nuclei interaction (from there the ; sign). Now we go back to 
the total problem 3.8 and use the expansion ^'(r, 91) = Ylk 
We should not enter in the specific details of calculus, just mention that the 
electronic pseudo-solutions are orthogonal, therefore one can integrate over 
the electronic coordinates the resulting equation and obtain, formally: 


[Tn + Ek- E]xkm = Y,Fk,iXi{^) (3-17) 

i 

= - E (3.18) 

Now, reminding that Mj ^ m and in general \Ek — Ei\ ^ 0 we could 
safely assume that the right hand side terms are small compared to the left 
side ones, therefore, we can write an eigenvalue problem also for nuclei. But 
this was the whole purpose from the start: to separate the total Schrodinger 
problem in two eigenvalue problems for electrons and nuclei, coupled only 
parametrically trough the interaction potentials. 

The BO approx, has common points with the adiabatic theorem |60j 
and some extended discussion can be done on the so called potential energy 
surfaces, given by the solution of the clamped Schroedinger eq. 

In clusters, the Born-Oppenheimer approximation holds almost indef¬ 
initely, even though it might be questionable in the ultra intense laser 
regimes. But in general system the problem is not that fortunate and one 
could refer to specific articles that treats [H], |62] the molecular problem 
beyond BO approximation. 

3.2.2 Dynamics of nuclei 

Now that we have seen how the full nuclei-electrons problem in a cluster can 
be separated in a problem of electrons and one for nuclei, let us go further 
and see what can be done more to simplify the treatment. In quantum 
chemistry, it is known by the name of Molecular dynamics the study of ionic 
geometry (in ground state or dynamics) with methods which involve a more 
or less detailed level of quantum effects. These studies coupled intrinsically 
the motion of electrons in the problem. In contrast with that, the present 
subsection discusses just the motion of nuclei(ions) in a cluster. 





To avoid any confusion, the BO approximation does not state that the 
nuclei and electrons are independent, in the sense that one could solve their 
problem separately, but rather that the problem is separable in a math¬ 
ematical sense of nuclear and electronic variables. As we shall see soon, 
their dynamics is strongly correlated with the interaction potential which is 
dependent on both systems. 


Starting from 3.17 Schroedinger problem for nuclei in the frame of BO 
approximation, we extend it to time dependent processes: 




+ {Vnn + Vne)x = i^dtX 


(3.19) 


We write some kind of Madelung transform [63] for the total wave func¬ 


tion of nuclei x(51, t) = A{d\,t)exp{iS{yi,t)/h) and separate the 3.19 equa¬ 
tion in real and imaginary parts, obtaining a continuity and a phase equa¬ 
tion: 


I 


dtA{9\, tf + Y t)/2Mi) 

I 

2Mi ^ + ynn + 2^ A(91, t) 


0(3.20) 

0(3.21) 

(3.22) 


Now, since it is known that the nuclei have diameters of fm order in 
comparison with our range of interest which is of A order, we can approxi¬ 
mate them in the classical limit of point like particles and go to the classical 
limit ^ —7- 0. Furthermore, let us denote the total potential interaction en¬ 
ergy between nuclei and nuclei-electrons with 17(91). Also, without entering 
in details, the remaining part of the equation for phase S (apart from the 
time derivative) can be written as a Hamilton function in corespondence 
with the Hamilton-Jacobi formulation of classical mechanics, and denoting 
the momentum of the I — th nuclei with Pj = V/5(91, t) we obtain the 
classical equations of motion for nuclei: 


Pf = P, (3.23) 

(3.24) 

As a conclusion of this particular subsection, we have obtained (under 
the assumption that we know the interaction energy of nuclei) the classical 
Newtonian like equation of motion. Of course, the potential energy V (91, t) 
depends not only on the Coulomb energy between nuclei positions but also 


dPi 

dt 








incorporates the electron-nuclei interaction which can be a quite disturbing 
term, depending on the method used to describe it. 

In the dynamical regime, supposing that we know the initial conditions 
for nuclei and electrons, we have only to solve the equations 3.23 with some 
associated numerical methods, again, presuming that we know at any time 
the Ven part of V. The problem of initial condition is rather hard in the 
sense that an usual type of cluster simulation starts with the system in its 
ground state and acts at t = 0 with an external laser field. Therefore, it is 
imperative to have full knowledge about the ground state configuration. 

The stationary cases of the equations of motion does not tell us anything 
new, just that the individual momentum nuclei must be null and remain that 
way, which means again that the configuration together with the electronic 
configuration must give the global minimum of V which is a functional. 

Being an optimization problem which is by default highly non-linear, we 
should think from the start to tackle it with an iterative method. Beside 
the iterative aspect, there are many numerical methods to solve optimization 
problems in the world of mathematicians. Nonetheless, not many of them 
apply to a molecular problem, in particular through this thesis, Monte Carlo 
simulated annealing methods have been used. Generally this method is 
slower than others, but in the problem of molecular optimization is necessary 
since the potential energy has some special features: the convexity is not 
known a priori, the number of variables is large and the number of local 
minima is also very large. Since we search for the global minimum, it is 
necessary to use a method which is able to tunnel through the local walls 
in the parameters space around a local minimum. More details about the 
numerics will be given in[^ 


3.3 Hartree Fock theory 

Historically, Hartree-Fock (HF) theory, appeared in |64j with the work of D. 
R. Hartree which assumed that the total wave function of a system of iden¬ 
tical particles can be written as a product of single particle wave functions. 
Obviously, this assumption is wrong since it neglects the antisymmetric char¬ 
acter of fermions therefore, the Pauli principle. His work has been refined 
later by Fock [65] and Slater [66| which took into account the antisymmet¬ 
ric feature of the total wave function. The equations which were obtained 
under the energy minimization principle were concerning the ground state 
of a system, but soon [55] Dirac proposed a time dependent extension of the 
theory. 

Later, new representations or formulations (second quantization, many- 
body theory, quantum-electrodynamics, path integral, quantum field theory) 
have been invented to deal with the quantum mechanics of many-body sys¬ 
tems and HF theory has been derived in many other ways, more different 



and rigorous than the original theory. 

To continue the logical path started with the QLiouville’s equations, we 
shall derive HF theory from density operator considerations. Recalling the 
BBGKY hierarchy |3.4t we take now the equation for the first order density 
operator 


Q'(f2) 


ihdtti = [FSfi] + Qi(f2) 

(3.25) 

i=k 

(3.26) 

{N - l)Tr2mxk - X2),t2]} 

(3.27) 


Now, there is a standard abstract form (named cluster expansion [67j) for 
the relation between two consecutive density matrices which for the first or¬ 
der, reads: f 2 (ri, r'^, r 2 , r^) = f i(ri, r'jf i(r 2 , r^) - f i(ri, r 2 )-^ 512 - 

We have used the position representation and the ^'12 is called the corre¬ 
lation operator. Now, the hrst step in the HF approximation is to neglect 
correlations with higher order density matrices, e.g. gi 2 = 0. This condition 
is well described in the quantum field theory by means of Feynmann dia¬ 
grams where neglecting this type of correlations is equivalent with a mean 
field theory, aspect which is essential in HF. 

If there were to be no spin, the first eq. from eqs. 3.25 could be simplified 
to ihdtTi = [H^ + [///,ri], where Uh is the mean field Hartree operator, 
Uh = T^’(<?^'i 2 r 2 ). But, taking into account the Spin-Statistic theorem [68] 


one can prove that the Hartree operator is defined by Uh = Tr((/)i 2 , f 2 ) 
where A± is the anti-symmetrization operator. Further let us write down 
the equation obeyed by the Fi and define the Fock operator: 


ihdtt, = [F, fi] (3.28) 

F = H^ + Uh (3.29) 

F{x, x') = {— -h n(x))5(x — x') -\- 6{x — x') [ - U'{x”, x")dx'' -r(x, a:()1.30) 

2m* J xi 2 xii 

Furthermore, the assumption that the second order density operator can 
be developed in a difference of product of first order density matrices is 
equivalent with the fact that Fi is idempotent. In turn, the idempotency 
implies that there is a natural set of eigenvalues jV'j) for Fi such that Fi = 

Yli (Note that idempotent means pi G {0, l},Vi G N). 

In the coordinate representation, this condition can be expressed through 
the fact that the total wave function is a Slater determinant and so, we 
have obtained through a factorization and an approximation, the original 




assumption used in the HF theory. With the above mentioned relations, the 
total energy in a HF system can be expressed only in terms of Fi : 


Ehf[^i] 


= / F(xi,X2)Fi(xi,X2) = y (-^+ n(xi))7i(xi,x'i)|,^^=,^ 
1 f f 7i(x2,X2)7i(xi,Xi) - 7i(x2,Xi)7i(xi,X2) 

dvreo yy |x2-xi| 


/dxi ^3.31) 
(ixi(i?^g5-32) 


Now, if one uses the minimization principle SEhfIT] = 0 with the con¬ 
strain rr[Fi] = N, obtains after some algebra, [F’, Fi] = 0, which is kind of a 
trivial information, since we knew that the ground state is a stationary state, 
thus, the Ihs from eq. 3.28 is 0. But the meaning of this goes further. The 
fact that the Fock operator and the first order density operator commute, 
means that they have common eigenfunctions, i.e. the {IV’i)}- Therefore, we 
could write down, using 3.31 the expression for energy in terms of orbitals 


i 

// ^ -rV'j(x2)V'j(x2)dxidX2 - 

47reo ^ J J X2 - xi 

// -rV'i(x2)V’*(x2)dxi(ix2 

X2 - xi 

Now, we use the constrains of = Sij with the energy mini¬ 
mization principle 6{Ehf — = 0 to obtain the stationary 

HF equations. Similarly, we can go to the time dependent regime where 
ihdtTi = [FH-i?,?!] (or, using in stead of 5E = 0, the principle of action 
extremization 5A = 0 where A = f dt{Y^j{ijjj\dt\'ijjj) — Ehf))- 


(3.33) 

(3.34) 

(3.35) 


EHFl'ipi) = ffilV’i) (3.36) 

Fhf\A) = iFdtltpi) (3.37) 

There are further details, that should be discussed regarding the N- 
representability, the orthogonalization method, but this kind of discussion 
could drive the present subsection to an unnecessary length. The main 
points to be retained from this are that (TD)HF theory is a mean field theory 
that neglects statistical correlations between first order density matrix and 
higher order matrices and it is representable through a single particle set of 
Schroedinger like equation. This equations contain a natural term that gives 
the Coulomb self interaction of the electron density and a supplementary 
non-local potential known as exchange. The latter one is a pure quantum 












effect arising from the anti-symmetry condition for the total wave-function 
(density matrix) and, as we shall see, is the main reason for which HF 
involves a greater computational cost than DFT (for example). 

It is worth mentioning that a true HF simulation should take into account 
the fermionic character of the orbitals beyond the anti-symmetric feature of 
the Slater determinant, including explicitly the spin. This is done working 
not with trivial wavefunctions but with spinors. 

As in any section of the present chapter, only the basic notions of the 
method are discussed. In practice, HF has a history of nearly 85 years in 
which has been used preponderantly in nuclear physics. Since many systems 
have been found where the results were not accurate enough, extensions (the 
so called post HF theories) have been invented. The most straightforward 
one si the Configuration Interaction (Cl) method which relaxes the condition 
of a single Slater T determinant for the total wavefunction to an expansion 
in a basis of Slater determinants T = Further, one can use a type 

of perturbation theory to extend the zeroth order Fock operator (eq. 3.30) 
to a ’’perturbed correct” Hamiltonian and so arrive at the Mpller-Plesset 
(MP) Perturbation Theory. Further, the Coupled Cluster (CC) theory uses 
an exponential cluster operator to improve the results. And the list can 
go on. The main purpose of all this extensions is to capture the electron 
correlation which was discarded from the start in the usual HF, since this 
quantity can be quite important in various systems. 


3.4 Density Functional Theory 

In the past decades. Density Functional Theory became the master method 
in many body simulations in a wide range of systems. Extensive studies can 
be found in nuclear physics [69], atomic physics, molecular [ZOI, cluster m, 
quantum plasma m , quantum dot m. solid state, etc. 

As we shall see, there are serious similarities with the HF theory, but 
there are also some fundamental differences which makes it faster in numer¬ 
ical simulations. This is the main reason for which is chosen over HF. 

Historically, the first genuine DFT theory appeared in the work of Thomas 
and Fermi IZU soon after Schroedinger equation was derived. They basically 
assumed that the density of kinetic energy of a fermionic system can be ap¬ 
proximated with the only density dependent expression derived analytically 
from the Homogeneous Electron Gas (HEG) model. This approximation 
gives, from a statistical perspective, an equation of state in the thermody¬ 
namic limit of large number of particles N ^ oo. We shall discuss in detail 
the TE approximation in a future section. 

Some other extensions as the Bloch model which is just the time depen¬ 
dent versions of TE have been worked out in the next years, but essentially, 
it remained with the status of a simple, not reliable model, until 1964 when 



in a paper m was putted the idea of DFT on solid mathematical ground 
with the two theorems known nowadays as HK theorems. One year later, 
m designed a more practical way of implementing the DFT with a set of 
non-interacting particles obeying the KS equations. From that moment the 
road was free to extensions and applications. In the 90's have recreated 
the work from the 60's for time dependent phenomena and in coherence 
with the development of the computers, it became the tool for world wide 
scientists in a lot of domains. 

The central idea of DFT is that ’’everything” can be done just knowing 
the density of particles. In terms of Fi, the density can be expressed as 
the diagonal part : p{r) = ri(r, r). While the entire original construction 
was done on grounds of p, assuming that the external potential in which 
our system is placed is purely local u(r) we shall take into consideration 
the fact that there are systems in which the the potential can be non-local 
u(r,r'). The latter case has been treated by m- For this reason we shall 
pass through the HK and G’s theorems in parallel. 

3.4.1 Hohenberg Kohn theorems 

It is a logical (and mathematical) fact that having the N number of particles 
fixed and the external potential u(r) all the properties of the ground state 
(GS) are fully (uniquely) determined. In DFT this idea is somehow reversed, 
in the sense that the external potential is uniquely determined by the density. 

Now, to state the first HK theorem and the first Gilbert theorem : 

HK: Between the external potential u(r) and the density of particles p{r) 
there is a hijective mapping in the sense that the density determines the 
potential up to a trivial constant. 

G: The external potential determines uniquely the density matrixTi{r, r') 

To sketch the proof, let us presume by reduction at absurdum that there 
are two total wave functions IT) and |T') and their associated external 
potentials u(r), u'(r) give us the same Fi and consequently the same p. If 
this is true, the ground state energy can be written for the two cases as: 


Egs = (3.38) 

Kgs = (T'|.H'|T') (3.39) 

Now, since all the properties of the GS are defined by the densities, one 
can separate the energy functional in an universal unknown functional of 
density and the potential energy given by the external potential Egs = 
F[p\ + f p(r}v(r)dr or Egs = -^[hi] + / F(r, r')u(r, r')firdr' for non-local 
potentials. The Ritz variational principle says that the ground state energy 
is minimum, therefore: 


Egs - E'as = j iir(ir'ri(r,r')(w(r,r') - ^;'(r,r')) < 0 

E'qs-Egs = j drdr'Ti{r,r'){v{r,r')-v'{r,r')) <0 

But, adding this two relations we obtain Eqs — Eqs < — Eqs which 

is absurd therefore, the hypothesis of the theorems hold true, both for local 
and non-local external potentials. 

HK: The true ground state density p{r) minimizes the energy funetional. 

HK: The true ground state density 7 i(r,r') minimizes the energy fune¬ 
tional. 

Now, having the same external potential we have the p{r) and 7 i(r, r') 
which minimize the energy functional E[p/'yi] and gives us the ground state 
configuration. If there would be another p' and an associated 7 ^ that would 
minimize the energy, then by the variational principle, those minima would 
not be the ground state, therefore, the density matrix or density of particles 
that gives us the minima in energy are unique. 

The Gilbert theorems have been used, to generalize the original DFT to 
non-local potentials, which, as we shall see later, are very common in cluster 
physics. 

Now, having this theorems that tells us that the Ground state density 
operator is uniquely determined by the external potential and is unique for 
the ground state, one can go back to the energy functional E and apply 
the energy minimization principle with the constrain of constant number 
of particles (in the microcanonical ensemble) to obtain an Euler Lagrange 
equation. Indeed: 


S{E[-fi] - pN} = 0 (3.40) 

d{E[p] - pN} = 0 (3.41) 

Using the explicit expression for the energy and the number of particles 
N = J p(x.)dx = ff 7 i(x, x')dx'dx we obtain the following Euler Lagrange 
equations: 


dF[p] 

5p 


-I- u(r) = p 


—- + v{r,r) = p 

oji 


(3.42) 

(3.43) 


The value of the DET is that, in principle, if one would know the E[p] 
or E[ 7 i] functionals exactly, then it would be straightforward to solve the 
above equations and to find the densities. Unfortunately, for this functionals 



only approximative expressions are known and will be discussed later. One 
might, wrongfully, think that could express the total energy as in eq. 3.31 
from the HF. But let us remind that in there some specific factorization of 
the r 2 in terms of Fi has been worked out plus an approximation of zero 
correlations. Here is not the case since one of the purpose of DFT is to 
capture as much as possible from the correlation effects. 


3.4.2 Kohn-Sham method 

In 1965 Kohn & Sham have partially cured the problem of unknown F 
functional in DFT (for local external potentials) introducing a system of non¬ 
interacting particles. In essence, one could separate a kinetic energy term 
and an interacting energy term in F in such a way that F[p\ = T[p] + Vee[p]- 
Even trying to represent the density and the kinetic energy in terms of 
natural orbitals of 71 (which are not to be confunded with the ones from 
HF), one could at best write p = and T[p] = 

But in here, the summation goes over an infinite number of orbitals. 

The idea of H&K was to take a particular case of this representation, 
in such a way that pi = l,Vz < N and basically reframe the problem to a 
set of non-interacting particles. Still, the kinetic energy provided by this 
set of orbitals is not the true kinetic energy, nor the interacting energy Vee 
can be exactly reproduced. Therefore we introduce Ts[p] = the 

kinetic energy of the non-interacting KS orbitals and p = With 

this, one can formally rewrite the total energy as Fhk[p] = Ts[p] + T[p] — 
Ts[p] + Vee[p\ + J[p] - J[p] = Ts[p] + J[p] + Exc[p]- Where the J[p] is the 
Coulomb interacting term like in the HF theory. 

So what?, one might ask. Now the energy is even more complicated since 
you have the unknown Exc[p] plus a term which is at best representable by 
a system of hctitious particles. Well, now we will take into account the min¬ 
imization principle for the total energy, with the associated normalization 
constrains for KS orbitals and performing the minimization with respect to 
orbitals. By some simple functional derivatives, one obtains the KS equa¬ 
tions for a ground state system: 


[~ ^ F VKs\'^j = (3.44) 

VKS = Vext + VH + Vxc[p] (3.45) 

To decript the terms from the effective KS potential vks we see that 
the first is obviously the external one, which has to be specified for every 
system in part (in particular, for clusters on the ground state, is the potential 
created by the nuclei or ions). The second, is called Hartree potential in 
connection with the effective potential which can be found in HF theory and 
is subject to a Poisson equation (having Coulomb nature) V‘^vh = —47r/9(r). 




Of course, talking about non-interacting particles, the total density can be 
easily reconstructed as p{r) = The last term, which stands for 

exchange-correlation potential, can be written formally as derivative of the 
exchange-correlation energy Vxc = ^^xd^P-, but apart from this we do not 
know its specific form. 

Now it should be clear that as HF, DFT allows us to solve N single par¬ 
ticle equations in the mean field approximation from which any observable 
can be computed. The point in which the DFT becomes easier to be used 
than HF is exactly the term which is unknown. There has been a lot of 
search for good approximations on Vxc and usually, it has a pure density 
dependent form, more or less simple or local. Still, being density dependent 
the calculation of this potential is much more easy to perform than the effect 
of the exchange operator from HF, therefore, the whole method is easier to 
implement numerically. 

Regarding the approximation used for Vxc we do not enter in this subject, 
since it is only of practical use for cluster physics, but much of the physicists 
or chemists that work in the DFT branch are drawn in the search for better 
functional. I will just remind with appropriate references the well known 
approximations. First is the ideal approximation, LDA m in which xc 
potential is local in the total density p therefore, very easy to compute 
numerically (generically the exchange part is taken as for the ideal gas oc 
An extension of LDA is to consider also local functional for xc but 
which include further gradient corrections and this approximation goes by 
the name of GGA [80]. Going even further, one can use functionals which 
include laplacian of the density, MetaGGA or go the the non-local potentials 

m- 

Through the simulations done for this thesis, I have used the Gunnarson- 
Lundqvuist [82] due to its LDA form. Still, there are problems with such 
approximation that must be corrected by the so called Self-Interaction Cor¬ 
rection [83] . 

As a hnal remark, DFT should incorporate (like HF) the spin character 
of the orbitals. This can be including working with spinors instead of wave- 
functions, or simply assigning a spinorial a label to each —)• The 

effects is that in practice there are xc functionals that take into account 
separated densities for spin {p^, pfi) and this can give certain differences in 
the numerical results, in particular on the total energy of the system. 

3.4.3 Time Dependent DFT 

Now we have seen how the DFT for local and non-local external potential has 
been constructed historically, we should go further to the time dependent 
version of this theory since, by default, the interaction with strong laser 
pulses is dynamic phenomenon. The basic work in this aspect has been 
carried out by Runge & Ross in 1984 m and essentially follows the same 


logic as the HK original DFT derivation. For this reasons we will not enter 
in too much detail, just perform a short parallel. 

RR:The bijective mapping between r;ea;t(r, t) and p{r,t) (or j{r,t)) holds 
for time dependent systems in the sense that the eternal potential is deter¬ 
mined uniquely up to a function of time. 

RR: The true density p{r, t) minimizes the quantum action functional. 

The proof of the first theorem is more involved than in the stationary 
case, but the second one is applied in the same manner: a functional called 
action A[p\ is defined and separated in a single particle manner which in turn 
gives us a set of Schrodinger like equations in the time dependent regime: 

ihdti(j = (3.46) 

While formally, we have a Schrodinger equation and the effective KS 
potential is constructed in the same way as in the stationary case, the time 
dependent DFT problem is a bit more involved. The first heavy impediment 
is that in the dynamics, many of the functional approximations for Exc are 
not valid anymore, since they are derived under ground state considerentes. 
Therefore, the field of time dependent Vxc approximations is a bit more 
dry, or at least gives worse results. A true functional should have memory 
properties in the sense that should be a time integral over a non-local time 
kernel, which would make any computational analysis far too complicated. 
For this reason, the LDA/GGA functionals are usually transferred in the 
dynamic case, at least for gross properties. 


3.5 Phase space and Vlasov limit 

Until now we have a picture of how from the density matrix formulation of 
quantum mechanics using Liouville von-Neumann equation and the reduced 
density matrix approach one can derive the BBGKY hierarchy. Further 
this hierarchy can be truncated at the zero level neglecting correlation and 
obtaining the HF. 

Now, we should go further on another approach of parallel power with HF 
but contained in a different representation of quantum mechanics, namely 
the phase space representation. The pylons of this direction have been 
putted by Wigner in its seminal paper |84) on the quantum effects at ther¬ 
modynamic systems. From there, a lot of work has been performed and now, 
there is a solid literature and mathematical apparatus that can be invoked 
in this direction. We shall start with a basic introduction in the formalism 
of the quantum phase space and than the Vlasov equation will be derived 
as a semiclassical limit. 


3.5.1 Phase space representation and Wigner-Weyl trans¬ 
form 

The search for alternative representations of quantum (statistical) mechan¬ 
ics has been always present in physics due to some intrinsic hope that a 
different representation would give access to the same physical reality from 
a different mathematical perspective. This is the case of the phase space 
representation which holds many similarities with the classical phase space 
statistical mechanics. 

Soon after Schrodinger’s equation, Hermann Weyl found a way to map 
the classical functions from the PS to operators through a quantisation pro¬ 
cedure [85] called Weyl quantization. Conversely, in 1933 [?] found a way 
to map the quantum operators into classical-like functions in the PS. This 
bijection from quantum operators to classical functions is contained mathe¬ 
matically in the so called Wigner- Weyl transform. Let us consider a function 
/[<!)] and its associated operator $[/]. The transform reads : 


m] = 2 [ dye-^^Py/\q + y\m\q-y) 


(3.47) 

(3.48) 


The position-momentum operators from the Weyl transform obey the Lie 
Algebra with the associated commutation relation: [P, Q] = PQ — QP = 
—ihl. One could easily slip on the path of mathematical questions regarding 
the Weyl algebra [86], Moyal bracket m, etc. While being a theoretical 
thesis, still, it is not my purpose here to treat such subjects which have no 
practical value. 

Now, we have a quantum theory, from the BBGKY hierarchy and want 
to obtain a PS one, therefore, for us the Weyl transform is useless. We 


shall focus further on the Wigner representation which transform the 3.4 
equations in a PS ones. First of all, the s reduced density matrix P^, after 
being dragged through a coordinate representation and than subject to a 
Wigner transform, gives us the evolution equation in PS: 



(3.49) 

r, t) = ^Q[V,Ts+i] 

(3.50) 


The Q term is quite tedious to show, but we will particularize the form 
for /i: 





•A + iv.}A(p>-.) = i / 

Q = ^(a;)|[^^r|/i(r,pi,i) + J dr2dp2(/'(a:)|^;^^2;^^|/2/i2(r-,pi,P2,r2,t)(3.52) 

Now, for brevity and elegance in writing we define the Moyal bracket 
as a mathematical operation which is basically a composition law which for 
two functions f,g defined on a 2d Euclidean phase space like {f,p}: 

{{f,9}} = ^/(r,p)sfn(^(^r^p - ^p^r))9{r,r) 

Where the sense of the ■<— means that it acts on the left side of the 
expression, namely on /. If one, as in HF neglects the correlation and 
generally the coupling with higher order matrices, the so called Quantum 
Wigner equation is obtained in the compact form: 


dtf{r,p,t) = -{{f{r,p,t),H{r,p,t)}} 


(3.53) 


From now on, /(r, p, t) = /i(r, p,t). This elegant form is interesting 
from several perspectives. First of all, we have been able to recast the whole 
Hamiltonian operator inside the Moyal bracket. Second, the similarities 


between eq. 3.53 and the classical kinetic equation (which also can be written 


for Hamiltonian systems in terms of Poisson brackets) is astonishing. If one 
looks closely to the Moyal bracket and takes the h —>• 0 limit, the classical 
Poisson bracket is obtained and the motion of particles becomes a classical 
one. The reason for the differences between classical and quantum can 
be understood in terms of phase space in the sense that the commutator 
Ip, x] = ihl is preserved in PS as volume. 

Further differences can be understood from the properties of the Wigner 
distribution function /(r, p, f). First of all, it is a real function which has 
its norm (the integral over the PS) equal with trace of the density oper¬ 
ator, so equal with N the number of particles. The evolution equation is 
space-time symmetric and also Galilei covariant. From the definition of 
the Wigner transform, any physical quantity A can be obtain integrating 
the product between /(r, p, t) and the associated function for that quantity 
Air, p). In particular, the density is the integral of f over the moment space 
Jdpf{r,p,t) = (r|p|r) = p(r). 

Perhaps the most important property is related with the limits of /. 
It can be proven that due to the structure of the PS and the evolution 
equation, —2/h < f{r,p,t) < 2/h. As one can see, / is not bounded bellow 
to zero which means that in can have local negative values. In this respect 
doesn’t met the classical criteria of probability distribution function and the 
interpretations of it can be quite doubtful. 




Even though, knowing / (in the frame of HE approximation) is equivalent 
with knowing basically any quantity of interest for our system, the path 
towards such knowledge is not just nontrivial, but usually impossible. The 
challenge of solving |3.53| even numerically is not feasible, therefore, further 
simplifications must be used. 

3.5.2 Vlasov equation 

In the classical kinetic theory there is widely known the Boltzmann equation 
for the distribution function. This equation has its quantum correspondent 
but both have a fundamental flaw: the long range forces (as Coulomb force 
between electrons) are not described. The extension known from plasma 
physics is the Vlasov equation. The quantum correspondent can be derived 
from the Quantum Wigner equation. 

If one takes the Moyal bracket and expand it in powers of h and then 
retains only the first two terms, the following semi-classical equation it is 
obtained: 


{dt + ^Vr - VrVVp - 

m 24m'^ 


V?VV3)/(r,p,t) = 0 


(3.54) 


Now, we have neglected in the derivation of Q Wigner Eq. 3.53 the 
rhs. of the zero order BBGKY equation. From a physical view, that part 
contained correlations of particles and terms associated with the so called 
collisions. If one would want to include such effect this can be done in an 
ad-hoc manner with a somehow empirical term in the Vlasov Eq. which 
should take care of the Pauli blocking. The first approach on this part was 
done by [88]. Using this, they wrote the so called VUU equation: 


{dt + ^Vr - VrUVp)/(r, p, t) = I[f{r, p, t)] (3.55) 

J J dp2dp3dp4Wi23A[flf2h] (3.56) 

^1234 = + P2 - P3 - P4)'5(p? + pi - pi - pI) (3.57) 


As one can see, the logical construction of the collision term is to block 
the presence of two particles in the same Heisenberg volume element of the 
PS: {2'kK)^ consistent with the uncertainty principle. The Vlasov equation 
has the advantage (beyond numerical) of turning a positive defined /. In 
contrast, if one includes quantum terms like in 3.54 this property is lost. 

Now, even assuming that the semi-classical limit is a valid approxima¬ 
tion, one would want to include exchange correlation effects which can be 
essential in some systems. For this, there is another way to derive VUU-|-xc 







using DFT. Going back to the Section [3.4.2 we see that one inherent assump¬ 
tion of DFT was that Fi = Performing a Wigner transform 

on this density matrix and using the Kohn-Sham Eqs. 3.46[ one can arrive 
in the semi-classical limit at a very similar equation with 3.54, but with a 
different potential: 


{dt + — Vr - '^rVks^p)f{r, p,t) = 0 (3.58) 

m 

Now, exchange-correlations are included in a mean field manner in the 
VUU+XC equation simply by the presence of the Vxc potential for vks- 
There are other approaches that try to include the correlations effect by 
relaxing the initial factorization of the second order density matrix but we 
do not discuss them here. 

Regarding the set of quantum effects recovered in this VUU scheme, we 
obviously can point out the XC and the Pauli-type collisional term. Beside 
that, there is another statistical quantum effects which must be included 
from the initial conditions of the equation: namely the initial value of the 
/(r, p, 0). As we shall see later, this configuration is consistent with the 
Thomas Fermi approximation and it reflects the idempotency of the density 
matrix of order 1. 

Still, complicated, from the mean field perspective, the Vlasov equation 
is now feasible for numerical applications using one of the zoo of methods. 
From there we mention just the test particle method |89| and the PIC [UU] 
codes. Moreover, there are even analytical results which can be drawn from 
this kinetic theory, from which the most stringent is the Landau damping. 

The theory of Landau damping is old |91j and mathematically involved 
so there is no reason to present it in here, but an essential aspect should 
be noted: by linearizing the most simple variant of VLasov eq, without xc 
or collisional effects (only the electrostatic), one can obtain a dispersion 
relation which shows that any wave in a Vlasov system will damp itself, 
in principle because of the dispersions of the velocities in the momentum 
space. Other classical interpretations show wave-particle like interactions, 
etc.. But it remains crucial (and it will be used as a numerical advantage) 
that any excited quantum state should have during the dynamics this en¬ 
tropy preserving phenomena which directs it towards the ground state. 

Not to be confused with the H theorem that shows that any asymptotic 
solution of Boltzmann equation will be driven finally towards a Boltzmann 
distribution. This is a problem since allows numerical propagations up to 
ps and must be complemented by a high number of pseudo-particle in order 
to overcome this numerical thermalization |89j . 

Regarding the validity of VUU model, this is driven mainly by the lack of 
quantum effects. At high temperatures, specific to strong laser fields (above 
lO^iL) the nature of the interaction is basically Coulombian with collisions 
and the phase space picture becomes fully valid, in the classical sense. 









3.6 Hydrodynamic models 


Even though we have this three possibilities of finding more or less approx¬ 
imative the dynamics of a system of fermions (DFT, HE and VUU) we still 
can be in trouble. There are very large systems (as we shall see in the 
numerical result section]^ which are not feasible to be studied with single 
particle methods (HE, DET). On the other hand, even the VUU can impose 
seriously numerical difficulties for such system. Or sometimes one simply 
needs to be faster in computation at the cost of precision in calculations. 

For that reason, we search for further simplification. In classical physics, 
the simplihcation of the Vlasov (or any kinetic approach) is called fluid 
dynamics. Etymologically, the name was assigned because the equations 
obtained by moments of a kinetic equation have the same structure with 
the Navier-Stokes equations known from hydrodynamics. 

In quantum systems, the situation is the same. One can derive from 
VUU, DFT or HE the so called Quantum Hydrodynamic Model (QHM). 
This has been used extensively in the past years in [92l |93l |9H (9^ in appli¬ 
cations for metal clusters, fullerenes, semiconductors, etc. We shall describe 
shortly three ways of derivation it, which, even though give similar results, 
have some fundamental differences. 

Kinetic derivation 

We, as physicists must relate the theoretical work with some practical 
results and while for a mathematician the study of VUU like equations 
fullfills his soul, we have to obtain physical, realistic, results from it. For 
that reason let us see how some local observables can be obtained from 
distribution function. We list bellow, the density /?, the current J and the 
kinetic pressure tensor P: 


p{y) = 

j /(r,P,t)(ir 

(3.59) 

J(r,t) = J 

P/(r,pT)c^P 

(3.60) 

-P(rT) = [ /(r,P,Q(p - u) 

1 (g) (p — u)(ip 

(3.61) 


Where u = J/p is the global velocity field. Now, this are moments of 
the distribution function in the momentum space. The list can go on for¬ 
ever. Integrating 3.54 this way, we obtain an infinite hierarchy of equations 
which form the hydrodynamic model. We will write down only the first two 
equations because usually only this two are solved numerically and as the 
order of the quantity increases the complexity of the corresponding equation 
is higher; 



(3.62) 

(3.63) 


dtp + V J — 0 

dtJ + + P]+ pvv = 0 

Now, as we see, any n order equation, contains the n + 1 moment of 
/. It is like BBGKY hierarchy all over again. The solution is to find a 
similar truncation at some point in which the n + 1 moment is described 
by the previous n ones. Usually this is done taking the thermodynamic 
limit N —>■ oo and consider an equation of state for the pressure. The most 
common one is the polytropic (Thomas Fermi) P = P[p]l, where the scalar 
pressure P[p] oc p'^ but other approximations can be used too. 

This model has some issues: first of all, the Landau damping is no longer 
present since the features of the momentum space are lost. So linear waves 
miss the natural kinetic damping. Second, the pressure tensor is unknown. 
Third, even if the pressure tensor would be exaclty known in terms of density 
and current, where are the quantum effects? The truth is that deriving it 
from Vlasov, the QHM contains the quantum effects in the pressure term 
in some obscure manner by the statistical properties of /. A more detailed 
discussion on the derivation can be found in [96l |98] . 


DFT derivation 

Regarding the derivation from DFT, there is not too much to say. Al¬ 
ready we have a quantum hydrostatic model, simply from the Euler -Lagrange 
equation T 


9[p\ + VKs[p] = P 

But further, if one looks at the derivation of TDDFT, sees that a set of 
QHM like equation are obtained m as the continuity equation dtp+V J = 0 
and a momentum equation dtJ = P = {[H,j]). 

Already, the continuity equation is not a surprise anymore, and we are 
able to obtain it from any angle. The fundamental issue is with the current 
equation. By default, the m equations state just that the evolution of the 
current field in a system of N identical particle should be equal with the 
expectation value of the Hamiltonian-current commutator. Again comput¬ 
ing what you can from the commutator, the same formal system of QHM 
equations is at hand: 


dtp + S/ J = 0 
dtJ + + P[p]] + pVVks = 0 


(3.64) 

(3.65) 







Where the tensor P is again unknown. As for the Vlasov derived QHM, 
one can make the assumption that the tensor is isotropic and can be related 
with the Tg the kinetic energy functional by P = ts[p]l. 

Single particle derivation 

This third approach can be found in literature [99] and can be derived 
in the frame of DFT or HF since they both work with single particle or¬ 
bitals. It has the main advantage that reproduces pure quantum effects. 
Let us consider the set of HF/DFT solutions {V'i}- We use the Madelung 
transform which is nothing else than a polar form of an wave function 
V’fc = Pf! exp{iSk/h) in terms of single particle density and phase. Then 
we introduce it in the HF/DFT equation and separate the real and the 
imaginary parts resulting, in order to obtain a set of two equation known 
as Madelung equations and are the first attempt to construct a quantum 
hydrodynamic theory: 


dtSk + 


(V^ 

2m* 


+ V 


dtPk + '^jk = 0 

^2 y2pl/2 ^ ^ 

2m* p^P 


(3.66) 

(3.67) 


Now, from the definition of the current operator is easy to show that 
the single particle current = pkVSk- Also we can associate with this 
the single particle velocity field Uk = VS’fc. With these, one arrives at the 
density-current hydrodynamics for a single particle: 


dtjk + V(^^^ - ^pV ® Vlnpfc) + pkVV = 0 (3.68) 

Now we use the additivity conditions for density YlkPkPk = P for 
currents '^kPkjk = J plus the notation of total velocity u = Jjp and add 
the eq. |3.68|in the following form: 


(9t/9-hVJ = 0 (3.69) 

atJ +V(^^^^) + Vn + pfcVV = 0 (3.70) 

P 

f^2 

n = '^PkiPkiuk -u)® {uk -u) - -^pkS/In Pk) (3.71) 
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Again, we have obtained the equations of QHM. Now it is nice to observe 
that kinetic energy tensor 3.71 contains the Bohm like terms which are 
pure quantum. These type of effects were not present in the definition of 
pressure tensor derived under Vlasov equation since the li —0 limit has been 
taken. Beside reproducing some quantum effects related with the gradients 











of density, this form of QHM has the same issues of not knowing explicitly 
the tensor. Further ways to approximate it will be discussed in the next 
section. 


3.6.1 Kinetic energy functional/quantum pressure 


Until now we have stated about DFT that its conceptual power lies in the 


fact that is exact in principle, written in Euler-Lagrange equations as 3.42 


Also, any of the forms of the QHM model derived above are also exact in 
principle since a QHM is nothing else than the time dependent version of the 
abstract DFT. Nonetheless, one still doesn’t know the exact xc potential, nor 
the kinetic energy functional/kinetic energy tensor. Trying to approximate 
these terms (without going on the path of KS equations) we should find ways 
to express Ts[p\ or more general the n[/9] functional. To gain more insight 
about how to make such approximations, let us relate the above mentioned 
quantities to the KS orbitals, or to the hydrodynamic variables {pk-,Sk}'- 


= = ( 3 - 72 ) 

'' k '' k 

^ ) (3.73) 
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n(Di) = '^PkiPkiuk - n) (g) {Uk -u)- i^^PkS^ X Vlnpfc) (3.74) 

k 

In literature there are a quite various set of approaches to approximate 
this terms. First of all, historically speaking is the Thomas Fermi approx¬ 
imation [] which is based on the assumption that in the ground state the 
orbitals behave like free non interacting particles, therefore obey the Fermi- 
Dirac statistics. From these assumptions, one can easily prove that the local 
density and the local density of kinetic energy in the ground state can be 
related by: 


tTF[p] = Fop^^^ (3.75) 

Going further, Von Wieszacker proposed a correction to TF approxima¬ 
tion which takes into account the gradient corrections to the functional: 

V2pV2 

tTF+w[p] = tTF[p] H- 'pi/^ (3.76) 

Further there is an extension of all this, with the Bloch density matrix 
in the frame of Wigner-Kirkwood expansion which gives us the so called 
gradient expansion functional and it can be written up to the fourth order 
in h as: 
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With kq = 3 / 5 ( 37 r^)^/^, K 2 = 1/36, K 4 = (6480(37r^)^/^) ^ for 311 sys¬ 


tems. 


Nowadays, the most employed orbital free method is the one that uses 
an exchange-correlation potential with the TF functional corrected by a 
coefficient dependent Weiszacker term. Pauli terms, linear response based 
functional, etc. are also useful methods to derive approximations for the 
kinetic energy of a KS system. 

For the dynamic functional there are some attempts to derive it but the 
results are not impressive. Essentially, the TDDFT theorems assure us that 
there is indeed such a tensor only density dependent, fact which allows one 
to write Ilfp]. As described in the previous section, many times the only 
thing that can be done is to consider that some functional for the kinetic 
energy in the ground state holds also in non-equilibrium situations (which is 
not true, since by definition an equation of state is designed for equilibrium) 
so the mathematical approximation is done as IT = ts{p)\. This is connected 
also with the kinetic theory in the sense that is implied that the tensor is 
diagonal and there are some symmetries in the phase space. In particular 
the geometrical structure of the distribution function is only translated in 
the momentum space. 

Some approximation can be done on the Bohm like part of the 3.71| 
tensor supposing that all the orbitals have approximately the same spatial 
amplitude. This approximation has been used in |93j . but it has no solid 
justification. The only systems in which it could work are the systems with 
quasi free particles. Imposing it heuristically one can use it in the following 
sense: 



(3.78) 


k 


In general, even though the single particle velocities are irotational V x 
VS* = 0, the total velocity field can have rotational feature V x tt 7 ^ 0. But 
some anzatz can be done that the rotational feature of the total velocity 
held is negligible, therefore u ^ VS. This assumption is not unrealistic since 
the natural excitation in a cluster has dipolar character and the rotational 
helds are weak. Moreover, in [95] there has been performed and analysis 
of the kinetic solution during the dynamics and it has been shown that in 
principle, the Fermi sphere from the momentum space is always shifted with 










an irrotational field. Further, using the approximation of diagonal pressure 
tensor IlcZas ~ pf{p)i the QHM equation are to be written as: 


dtS + 


2m* 


dtp + VipVS) = 0 


+ fip) + ^ + 


2m* 


= 0 


(3.79) 

(3.80) 

(3.81) 


As a final step, the pseudo total wave function ‘h = p^/'^exp{iSjh) can 
be defined as an inverse Madelung transform and the 3.79| eqs. can be 
embedded in a Schroedinger like equation: 


^2 

ihdt<^ = [ -+ w]^ (3.82) 

2m* 

w = V + f{p) (3.83) 

This approximation has been studied in connection with semiconductors, 
metal clusters [95], fullerenes [inOj . thin metal foils, with good gross results. 
For a more detailed discussion about the Orbital-Free DFT see m- 

3.7 Linearised approaches 

We, as civilization, have constructed the mathematical theories in patho¬ 
logical way, based on a linear thinking. Unfortunately, the nature proved 
to be mainly described by non-linear phenomena and many times we have 
to appeal to the last resort measure: the numerical tools. Still, wanting to 
capture weak effects or normal modes in non-linear theories, one often does a 
linearisation of the equation around some known solution to get insight into 
the problem. The many-body theory makes no exception and every theory 
from the one discussed until now have been dragged trough the method of 
linearisation. 

While from an operatorial procedure we could go to things like Dyson 
series, or to Feynmann diagrams in the many-body theory, in the present 
section I shall avoid such things, keeping a more practical view on how DFT, 
HF and Vlasov equations can be linearised into other tools for weak laser 
cluster interaction. 

3.7.1 LRDFT 

In quantum mechanics, the most basic linearization procedure goes by the 
name of perturbation theory and is easily understandable in the Dirac repre¬ 
sentation. For many-body systems one can go to more abstract perturbation 
series like Dyson series. In principle the linearization of DFT equations can 












be derived from many angles. In the present work, I shall adopt a derivation 
based on the fundamentals of the linear response theory. 

In the frame of DFT one can speak about the single particle response 
function xo(r, r^w). On the other hand, the true linear response function 
x(r, r'jw) can be related (5p(r,a;) is the induced density variation in the 
energetic space) to the external excitation by: 


5p{Y,uj) = J dr'xo{r,r',uj)VKs{r') 
6p{r,uj)= / dr'x{r,r',uj)Vextir',u}) 


(3.84) 

(3.85) 


One can introduce in 3.84 the expression for Vks and derive a self- 
consistent equation for the density in the spirit of Dyson series: 


5p{r,uj) = j dr xoir,r', ui)[Vext{r,oj) + j 

, , , ,V’i(r)V';)(r)V'/.(r')V'r(r') , ^/’r(r)V’M(r)V';)(r')V'i(r'),„ 

Xo(r,r,a;)= > -^^--h- -f -——— 

Si- Emu - nw -ip Si- Emu + nw + ip 

i£occ,^£unocc 

This equation must be solved iteratively for the density usually on a 
real spatial grid. A description and more details on the method can be 
found in [102] . Once you have the density in the frequency domain, one can 
easily compute quantities as the dipolar response to investigate the optical 
response of the cluster for example in the linear regime. 

Finally one can recover a formal solution for the response function which 
will be useful to compare with the one from Vlasov theory: 


X(r,r',a;) = [1 - xo(r,r',a;)^^^^^|^] ^xo(r,r',a;) (3.88) 

3.7.2 Random Phase Approximation 

The Random Phase Approximation has some controversal history. Origi¬ 
nally it has been derived for nuclear physics and condensed matter calcu¬ 
lation 110311104111051 llOBj . An older version of it can be found as Tamm- 
Dancoff |107j approximation. It can be derived from considerentes of HF 
theory, but there are some problems related with its interpretation in terms 
of Feynmann diagrams where it can be proven that RPA is in fact a conse¬ 
quence of the Ring approximation |108) . By abuse of language, now both 
are known as RPA. 

In the present subsection we will derive the RPA equation following [109] 
. We start from the HF theory in which the total wave function is thought 





















to be a Slater determinant j'l'o) in the ground state. Now, we act with some 
external perturbation which induces an rjG generator of the dynamics. By 
means of Thouless theorem m one can show that the time dependent state 
of our system can be written I'h) = exp{ir]G}\^ q), where G is hermitian, 
and moreover must be a first order particle-hole like operator, G ~ apah- 
In order to find the equations of motion for the system we apply the 
variation principle of quantum action A = f minimization 

5A = 0. We do this, inserting the expression of |'k) in A, expand it to 
second order in r/ and minimize. Denoting {dpdh, dh^p} = ^la, the equation 
of motion is obtained: 

('I/o|[ata,i9iG]|'I/o) = (d/o|[ata, [i^,G]]|^o) 

The approximation G = + C~exp{iuit)] is taken in 

order to find normal modes in a fermionic system. Here is the operator 
that generates/annihilates the n — th eigenmode. Now, introducing these 
into the eqs. of motion and separating the linear independent oscillating 
factors exp{±iujt), we get two eigenvalue problems for the n — th eigenmode: 


w„(To|[ata,C't]|To) = (To|[ata, [iJ, C't]]|To) (3.89) 

-w„(To|[ata,C„]|To) = (To|[ata, [iJ, (7„]]|To) (3.90) 

Even more, to arrive at the well known matriceal RPA , we expand 
Gn = '^phi^^ph^p^h. — y^^^\dp) and obtain the standard eigenvalue problem 

( \ _ f ^ ^ \ f \ 

y^n) J -B* -A* ) V ) 

Where Aij^p^ = {ip\V\jiy) and = {i^p\V\ji). Solving the RPA 

eigenvalue problem poses a serious numerical difficulty in the computation 
of the A&iB coefficients. A discussion on this matter will be given in the 
Chapter 

3.7.3 Linearized Vlasov 

RPA/LRDFT are methods derived under the single particle picture of the 
system of N particles with the purpose of determining the linear dynamics 
in terms of eigenmodes. If one looks at a optical spectra constructed from 
RPA calculations, it will see a set of lines positioned at the resulting exiting 
energies. And this is natural, since in the search for eigenmodes we have 
assumed pure oscillatory behaviour. Sometimes, as using the Xs density 
response function, instead of delta functions, the resonances are fitted with 
some narrow lorentzians, given a phenomenological damping imposed by a 
fake imaginary part of the energy. 



Even though essentially a poorer approximation, Vlasov method is able 
to reproduce some damping effects even in the linear regime in the absence 
of any collisional integral. The phenomenon is known as Landau damping 
and its proof can be quite mathematically involved. The derivation will 
be just sketched further. We start from the Vlasov eq. 3.55[ neglect the 
collision integral, and assume that in the stationary state the solution is 
fo{r,p,t). We take an initial perturbation f{r,p,t) = fo{'f',p,t) + r]fi{r,p,t) 
from where the linearized Vlasov eq becomes: 


idt+v)h-^ I ( 3 . 91 ) 

Where Vi (r, t) is only taken to be only the hartree potential created 
by the perturbation density pi{r,t) = / dpfi{r,p,t), so V^Vi = —d-Trpi. 
Now we use the Laplace-Fourier transform of this equation and after some 
algebra the relation between the induced density and the effective field can 
be written as pi{k,uj) = where it has been introduced the 

density response function (similar with the same quantity defined in the 
frame of DFT). 


n(A:, to) = h 


dp 


Mp) - Mp + hk) 


(27rfi)3 _|_ piI2m — {p + hk)‘^/2m + i6 


(3.92) 


The relation 3.92 is known as Lindhard polarization function and it is 


equivalent with the RPA polarization function. When one does the classical 
limit of this eq. (the long wavelength limit /c —)• 0) obtains: 


n(A:, to) 



j pdfoM 

dv k^§^ 

(27r/i)3 ^ 


(3.93) 


The apparent problem with this functions is its continuity in the complex 
plane. Landau solve this issue and the interesting result is what now is 
known as Landau damping. Just to sketch the idea, one should consider the 
dielectric function 11 defined in the complex plane and study its behaviour 
when the energy (in particular 5) passes through the real axes. Writing : 


n(A:, Lu,d) = < 


V 


U^{k,uj,5),5 < 0 

duj'Il(k,u') ^ 

-;-Z7rn(A:, ui), J = 0 

j 2 tt u — u' 

IT^(A:, u, 6) — 27rill{k, uj,5),5 > 0 


(3.94) 


Where n'^(A;,ca, J) and Ii^{k^oj,6) are the usual retarded and advanced 
functions. 














3.8 Classical models 


When one goes to ultra intense laser fields, the quantum effects start to fade 
away. This should allow in principle the use of semi-classical or even classical 
method for the study of dynamics. On the other hand, there is a numerical 
problem given the fact that pure quantum methods (as DFT is) require the 
propagation of continuous fields which in turn require refined grids to be 
represented (also a good representation of the continuum spectrum). The 
electromagnetic problem itself is a large scale (compared to the electronic 
scale) problem therefore, the numerical refinement becomes a too hard task 
for a computer simulation. A last point is the fact that DFT/HF are theories 
designed for pure states and can be taken into account in the statistical cases 
using impure states (combination of non-integer occupation numbers). This 
means that at non zero temperature they describe canonical ensembles which 
by default can not resolve the microfluctuations arising in the strong field 
and which might be important for the dynamics. 

For all these reasons, the cluster world entered in the domain of very 
strong laser fields with classical methods. In this sense are two approaches: 
the molecular dynamics which basically denies any quantum effects and 
solves a classical point-like problem both for electrons and ions and the 
nano-plasma model which has been introduce as a schematic model to in¬ 
vestigate the macroscopic quantities of a cluster during its interaction with 
the external field. 

3.8.1 Molecular dynamics 

In the present, it is almost unanimous accepted that the single solution 
to solve the problem of cluster dynamics in very strong laser fields and 
consequently to represent the atomic ionization is the Molecular Dynamics 
techniques. Various groups have been used this approach and related works 
can be found in [111111121111311114L IHUl 11151111611117j etc. 

The essence of the theory is to consider that due to large amplitude 
of the electric field from the laser pulse, the electrons from the core levels 
of atoms are authomatically driven in the continuum spectrum, therefore, 
they become equivalent with the valence electrons, i.e. free particles. From 
this moment on, the MD methods should work taking into account only the 
Coulomb interaction in the electron-ion plasma formed. The interaction is 
cut in the near region since it can give unphysical electron-ion recombina¬ 
tions: 
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(3.95) 


The equations of motion for electrons and nuclei are solved as classical 
Newton equations. Still, the inner ionizations pose a difficult problem since 









require large energies and very short time scales. For this reasons, they 
are not treated explicitly but taken into account using probabilistic rate 
equations for the rate ionizations. 

For clusters which contain a few hundreds of electrons and nuclei, the 
MD problem is moderate from numerical point of view. When you pass to 
large clusters with thousands of atoms the problem is not directly tractable 
any more in the sense that the interaction becomes very expensive if it is 
evaluated with a direct method (scales with n^). In turn one has to imple¬ 
ment tree-based methods [] or to compute continuously the electromagnetic 
field using Particle-In-Cell techniques. 

3.8.2 Nano plasma model 

The last level in the chain of approximation is the nano-plasma model de¬ 
veloped in |114j . Its main assumption is that the rapid ionization creates a 
quasihomogeneous plasma for large clusters at the level at nm. This spatial 
scale requirement is given by the fact that the Debye length — 

5A should be less smaller than the spatial extension of the cluster. In turn, 
this motivates the treatment of global quantities and not microscopic ones. 

The main quantities to be investigated are: Nj the number of ions with 
the charge j, Ne the number of inner-ionized electrons, Eint the internal 
energy of the cluster and R its radius. These must be complemented by the 
Wj, the ionization rate for ions in charge state j, Q the total charge of the 
cluster. Pc = j (Svri?^) the coulomb pressure and Ph = UekbTe the hot 

electron gas pressure. With all these, the rate equations for nano-plasma 
model can be resumed at: 



(3.96) 


(3.97) 


dR _ Pc + Ph 5 
dt niiTii R 


(3.98) 


The total ionization rate can be decomposed in tunnelling, laser and 
thermal rate: 








Wj = 

_ f32mjj^ 

^ ~ V 

T4^j“^ = ne{crjv)ias 

W^ = ne{crjv)th 

Of course, = ‘iN^/{^ ttR^) is the concentration of electrons, Ip{j) is 
the ionization potential for the charge state j and Cn^i^m is a factor which 
depends on the quantum numbers describing the electrons involved in ion¬ 
ization. The term involving the variation of total charge can be further 
decomposed in: 
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(3.104) 

(3.105) 

(3.106) 
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For the internal energy and its energetic coupling with the external laser 
field, we can write: 


dEint _ p 2Eint dR dNj f3 IQS') 

dt “ R dt dt ^ ^ 

j 

Pabs = -^£intIm[£{0Jias)] (3.109) 

Where €(cj) = 1 — -|- zjzo;), the Mie frequency tCp = jme 

and the Landau damping term v = Av / R. 

Even though heavily simplified, the nano-plasma model captures almost 
all the important effects that take place during the dynamics. Therefore, it 
should not be used as a state of art/precision model, but as a starting point 
for the gross properties and trends to be expected. 



3.9 Final words about electron-nuclei interaction 


3.9.1 Pseudopotentials 


All the methods descried above that should be used for electrons assumed 
that the electron-nuclei interaction it is contained in the external potential. 
Also, in Section 3.2.2[ the nuclei-nuclei interaction has been described gener- 
ically with an effective potential. Now is the time to give more details about 
those and to present further simplifications that can be used. 

Naturally, the first assumption (the correct assumption) is that since the 
nuclei are considered classical, their interaction can be modelled with elec¬ 
tromagnetic forces, i.e. the Lorentz force:Fen = <lE[p,j] + v x B[p,j]. The 
electric and the magnetic field are subject to Maxwell’s equations and de¬ 
pend on the total charge density and total current density in the system. So, 
during the classical dynamics of ions, the Maxwell’s equation must be solved 
to compute the electromagnetic field created by the electrons. Fortunately, 
in clusters, due to their form and structure, the total angular momentum 
and the magnetic self consistent field are very small, and in this approxima¬ 
tion, the only relevant aspect of the interaction is the electric, Coulombian 
one. This can be obtained as solution of the Poisson equation. 

The nuclei-nuclei interaction, must be treated as classical two-body elec¬ 
trostatic interaction and the force acting on the I — th nuclei due to the 
other N — 1 nuclei is : 


j^i 


Ri — Rj 
\Ri-Rj\^ 


Only in the hydrodynamic models that treat the nuclei as a continuum 
distribution [] this force is also subject to Poisson eq. Conversely, the exter¬ 
nal potential created by the nuclei on the electrons is calculated as : 


^.nuclei 

^Ext 
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r - Ri\ 


But it happens that in the atoms from a cluster the core electrons from 
the filled shells suffer almost no dynamics if the external laser field has 
moderate intensities. To gain some insight in why is so, let us just look at 
the energetic levels in the Na and C atoms which are around 40eP and 2QeV 
for the 2p core electrons. Therefore, it would be a major simplification to 
get read of this core electrons considering them fixed with the nuclei. This 
gives us access to a treatment of ions instead of nuclei. 

But now, there is a question of how to treat the remaining electron-ion 
interaction. At first site, if the core electrons are considered bound in the 
point like region of nuclei, it should be straight forward to simply modify the 
q charge of the an ion. But keeping in mind that even bound on small spatial 
regions, the core electrons still have a tail of density at larger distances due 







to their wave function, it becomes clear that the screening problem should 
be carefully treated. Anyhow, the resulting mean field potential created by 
nuclei-core electrons complex is called pseudopotential. 

The first simplified models that treat this problems are solid core model 
|118j that uses the anzatz of solid sphere of charge for the core elctrons. A 
more refined model uses a density of core e in terms of superposition of two 
gaussians charge distributions. This gives us a pseudopotential of the form: 


Vpsp(r) = vpsP (|r-Ri|) (3.110) 

/ 

vpsp{r) = Y, Cj I (3.111) 

j 

There are other extensions to generalize this simplified approach. Gen¬ 
erally speaking, a local pseudopotential as the one above, is good only for 
alkali metals when you need to treat only the valence which have zero angu¬ 
lar momentum. For other types of elements, a pseudopotential must reflect 
the angular character of the wave function. The practical way of describing 
the pseudopotential in a non-local form: 


Vpspir) = Viocir) + J] VLir)UL (3.112) 

L 

^l = Y [ d^^'YLM{^)YLM{^') (3.113) 

M d 

Further ’’ultra-soft” pseudopotential can be constructed by the use of 
some parametrized solutions of KS equations for the valence electrons in¬ 
cluded in the general formula of : 


^ |V’c)(£i; - ec)(V’c| (3.114) 
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3.9.2 Jellium model 

Even further simplification can be worked out for the pseudopotential. In 
general, a PsP breaks any symmetry and imposes a full 3D treatment which 
can be quite expensive. On the other hand, there is a lot of experimental 
knowledge |119] about the almost symmetric (spherically, or axial) clusters, 
therefore, an approximation should be used in this sense to impose the 
symmetry in terms of the PSP. The method is called ’’symmetry averaged 









Almost straightforward we define the spherically averaged pseudo-potential- 
SAPS (for spherical symmetric clusters) and the cylindrically averaged pseudo- 
potentials-CAPS (for azimuthal symmetric clusters): 

U^irir) = j VPsPi^en - Ri) (3.115) 

d I 

Uback^{r,z)= f d4>^Vpsp{{rcos{(f)),rsin{(j)),z} - Ri) (3.116) 

d I 

Now, as long as the effective pseudopotential is continuous and second- 
order derivable, it can be linked to a pseudo-density of charge by the trivial 
V‘^Uhack oc Pps- This idea is not just a mathematical thing, but has a good 
physical sense due to the fact that the PsP is an artefact of the complicated 
electric interaction and screening of the nuclei-core electrons complex. This 
charge can be also averaged in the same sense as the effective resulting 
potential to give access to so called jellium model. The model was used 
for the first time in the condensed matter physics (for metals) where it was 
considered that the electronic density and the ion density cancel eachother 
in an exact way. 

For clusters, metal clusters, it is customary to use the jellium model as 
some distribution of positive charge pj^i from which the external potential 
for electrons can be computed as V'^Vbackir) = —ATrpjei- In a more general 
sense, it can be putted in relation with a two particle PsP interaction in the 
form: 


Vk.. = /<ir'Vp.p(r 


For spherical, metal clusters a lot of success has been obtained with the 
step jellium model [] in which : 


Pjel PO®i.\Rjel{r)—r\) (3.117) 

Rjel = -Ro(l + ^ aZmTjm(II)) (3.118) 

Im 

Where, obvously, Fjm are spherical harmonics and aim = oi*i-m- 
order to improve the results including diffusion)] effects at the surface of the 
cluster, soft jellium models)] have been proposed: 

Pjel = /0o)l + —-)]~^ (3.119) 

CTjel 

The normalization condition must be satisfied by the jellium density: 
f drpjei{r) = Nw — q. Where q is the ionization of the cluster. 



3.9.3 Atomic forces 


The pseudopotential method is a good method to take into account the 
screening effects when neglecting core electrons and usually gives good re¬ 
sults regarding the structure of a cluster. Still, one might be less interested 
in the fine details of electronic orbitals or simply wants to work in intense 
regimes where the details are not important. For such cases, empirical in¬ 
teractions can be introduced to model ion-ion repulsion/attractions. 

Generically, this is the case of rare gase clusters. If such an inter-atomic 
force is considered, then it can be use to obtain both the numerical optimized 
geometry and in molecular dynamics simulations also. To mention just two 
such examples, we write down the Lennard-Jones and Gupta potentials: 


VlAt) = 4e[{^r - i^f] 

(3.120) 

Uciri) = 

(3.121) 
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Vij = exp[-p{^ - 1)] 

ro 

(3.122) 

Pi = '^exp[-2q(-^ - 1)] 

^ ro 

(3.123) 



Chapter 4 

Numerical results 


Now that the main phenomena specific to laser cluster interaction have been 
presented in the second chapter and the main theoretical tools for study in 
the third, it is time to present the numerical methods that have to be used 
in connection with these theories to reproduce experimental data. 

It is my intention in this chapter to present personal results from as 
many different theories and regimes as possible. While almost each one of 
them, has been done long ago and can be found in different studies (to be 
cited), this chapter is a reflection of the personal work done in this direction. 

As structure, I have chosen to present the results divided in regimes and 
quantities of interest and then, for every regime many methods have been 
applied. Where possible, the simulations have been performed on each of 
the four cases presented in the introduction: Na, C , Ar and Xe clusters. 

All the simulations have been performed on different calculus machines 
with computing powers comparable with a standard personal computer. De¬ 
pending on the type of simulation to be done and on the numerical methods 
to be applied, the FORTRAN programming language or the WOLFRAM 
MATHEMATICA system of computation were employed. As a general 
thumb rule, for methods that implied test particle methods, Fortran codes 
have been developed while for partial differential equations, Mathematica. 
All figures have been created with the astonishing features provided in the 
Wolfram system. 

While times of computation could be specified individually, in general, 
there has been a concious sloppiness in the implementation in the sense that 
there have been used the minimum numerical refinement (in spatial grids, 
time steps, etc.) in order to gain computational speed. As a wide picture, 
semi-classical calculations as Thomas Fermi for ground states or calculation 
of static quantities are usually subject of a few seconds of CPU time, while 
at the other pole, TD-LDA or Vlasov methods had even days of computing 
CPU time. 
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4.1 Numerical generalities 

Following chapter almost every theory discussed in there ends with one 
or more partial (integro)differential equations (PDFs). Sometimes there are 
even systems of such equations coupled through some macroscopic quantity 
as the mean field potential. The truth is that this type of PDF’s are the 
practical tool to investigate the realityl in any sub field of physics, therefore, 
there is a long history of numerical, approaches more or less general with 
different levels of accuracy. 

Trying to form a clear picture on the logic and the relationships between 
different numerical methods that must be used in the laser-cluster physics, 
it would be nice to understand the problem from a more general, abstract 
point of view. Following basic single particle quantum mechanics, one can 
easily understand a PDF from the algebraic perspective. More precisely, any 
function /(r, t) that we encounter in our eqs. is in fact an abstract vector 
I/) in a Hilbert space. Any operation that is done on this function in our 
eq. is nothing else than the action of an operator on the function associated 
with that Hilbert space. 

Being abstract objects, we don’t have any understanding on the vectors, 
unless we represent them in a natural space for us. The most common space 
is the coordinate space, described by the infinite continuum basis {|r)}. 
The same is done with the operators and eqs like D\f) = \u) become in the 
coordinate representation D{r)f{r) = u{r) where /(r) is to be understood 
as the scalar product between our state |/) and some element |r) from the 
basis. Operators like momentum, are also local in r space but they become 
differential operators (gradient, laplacian, etc.). The choice of coordinate 
representation is done usually because is a natural link with what we see 
in real life, but otherwise, is not necessary. Other representations as the 
momentum representation can be useful since it can diagonalize gradient or 
laplacian operators. 

In general, having the set called orthonormal basis |ej) in which we make 
our representation of the reality, one can speak about a ’’box of interest” 
B. This is the subset |ep from the entire basis on which we expect our 
quantities to have significant values (also it is motivated by the fact that 
we cannot work in practice with infinite number of elements). In terms of 
|r) basis, the box can be visualized like an actual box, cube for example, 
in 3D. Also, we can speak about the boundary of the box dB which is a 
subset of B with some specific property. For example: if \e*) are vectors 
associated with angular variables (spherical harmonics) the boundary of the 
box is represented by the vector with the higher and the lower angular 
momentum. Another intuitive example is the actual box in the space 
which can be the surface of a sphere, or the surface of a cube, etc. depending 
on the actual coordinate system that has been employed. 


4.1.1 Numerical representations 

Now, the most natural numerical method, namely the Finite Difference 
Method (FDM) can be understood as a representation of a PDF in the co¬ 
ordinate space taking just ’’some” position vectors from a box. In terms of 
functions is to be understood as an expansion in delta functions. When the 
chosen subset of elements has certain regularity, one can speak about ” struc¬ 
tured grids”. For example, a standard 3D grid for finite difference consists 
from a box B = [—L, L] x [-L,L] X [—L, L] where the chopped basis of di¬ 
mension is the set of {(xj, i/j, Zk)},i, j,k = l,n with Xi = —L+2{i — l)L/n 
and the state |/) becomes fully described by the set of coefficients In 

the frame of this representation, operators like gradient or laplacian can be 
approximated with stencil configurations: 

fi+l,j,k fi—l,j,k 
hx 

fi+l,j,k “1“ fi—l,j,k ‘^fi,j,k 
hx"^ 

A more general representation is used in HF or in DFT methods. Even 
though there are many approaches that use real space with Eulerian grids 
more sophisticated than the one described above, there is another way to 
solve equations like Kohn-Sham. Instead of choosing the {|r)} basis, one 
uses another finite basis \e-*). Usually this basis has nice properties with 
analytical behaviour regarding how different operators act on it and more 
important is chosen in such a way to meet the criteria of symmetry and 
asymptotic values for the specific system. I will enter later in more details 
about how this method is applied on KS eqs., but essentially, the final result 
is numerically the same with FDM. 

Just as a remainder, the Fourier/Laplace/Fourier-Laplace/Hermite/etc. 
transforms are all spectral methods that can and will be used since can make 
a lot more easier the life of anyone. 

All this type of spectral methods can be regarded as Eulerian methods 

because, even though may not be explicit, they have in the background 

a ’’fixed” grid representation. But there is another type of spectral-like 

method which is more related with the Lagrangian methods or particle based 

methods. Essentially, it uses a basis with some unfixed parameter which is 

free to change under time evolution. To be more intuitive, let us imagine 

that we have /(r) and we choose a basis which can be also r represented 

as a function W{r — Vj] h) where /i is a parameter to control how narrow is 

function IF. Based on the identity /(r) = f dr' f{r')5{r — r') and assuming 

that f IF(r; h)dr = 1 and lim IF(r; h) = d(r), one can approximate; 

hs-O 

f(r) Rs ^ rriiW (r - r*; h) 


dxxf{r) 


dxf{r) 



This approximation is heavily used in what is called Smoothed Particle 
Hydrodynamics (SPH) [120j where usually the kernel W(r; h) oc exp{—r‘^/h?) 
The same exponential is used in test particle methods (TPD) |121] employed 
in the Vlasov equation. The main advantage of this method is that the due 
to the structure of the PDE’s involved (and their linearity), one can show 
that the entire time evolution of the quantity f{r,t) can be reduced to the 
classical motion of the particles represented by the Wi. The single aspect is 
that this particles move in a mean field acceleration field created by them¬ 
selves accordingly with the type of interaction present in the system. 

4.1.2 Types of equations 

Poisson equation is the most pregnant appearance in all the theoretical 
methods discussed since it models one of the basic type of interaction, the 
Coulombian interaction which usually supersedes any other one (like xc). 
In essence, the equation can be written V^</> = Itt/?. Working under FDM, 
one can write (j) and p as vectors of dimension (where d is the dimen¬ 
sionality of the problem: ID, 2D or 3D) and the laplacian as a x 
matrix. Further, you get a matrix equation which should be solved: a) by 
inversion of the laplacian matrix, which even though is sparse is the worst 
approach ever; b) by iterative methods as Gauss-Seidel, Jacobi, Succesive- 
Over-Relaxation (SOR), acceptable methods (for reviews and extensive dis¬ 
cussions see [12211123] !. Nowadays, one can find out there an entire zoo of 
numerical libraries that solve the Poisson eq. so you don’t have to do it by 
yourself. Working in full 3D, it becomes very expensive to use FDM. There 
are other approaches like finite element [124] . etc[], but the one that I have 
used during the simulations is the Fourier Transform. Taking the Poisson 
equation and performing a FT, one obtains a natural diagonalized solution 
and actually the solution can be written formally as (j){r) = P~^[P[p{r)]/k‘^]. 
The boundary conditions (BC) are essential in the sense that if you do not 
impose the right BC, you may end up with incredibly wrong results. In 
general, the boundary of the box is far enough to consider that the entire 
charge is contained in the box, therefore, a multipole expansion is used to 
compute the BCs. 

Eigenvalue Problems appear natural in HF and DFT as stationary 
Schrodinger equations: Xi\fi) = H\fi). Beyond any particular representa¬ 
tion, let us consider that we take the spectral box as B = {\ej)} which is 
large, finite and not necessary orthonormal. We can expand more like an 
approximation and not as a consequence of the nuclear spectral theorem the 
eigenvectors as \fi) = X^jC^lej), where = {ej\fi). Now let us introduce 
the following notations:4Iij = {ei\H\ej), Sij = (ej|ej), Cj = {c^,i = l,Af} 
the latter being the vector of coefficients for the i — th eigenvector. Now the 
problem can be rewritten as a matrix eigenvalue problem XiSCi = HCi. In 
literature, this type of equations can be found as Roothaan-Hall equations 











Again, for this problem, many standard numerical methods can 
be used m- 

Time-dependent equations arise naturally when you go from ground 
state towards dynamic regimes and you have to solve equations like time 
dependent KS. Their general form can be expressed in abstract as idt\fj) = 
H\fj). Being basically a Cauchy problem on time, this eqs. have the formal 
solution: 


\fj{t)) = exp[-i/h j dTH{T)]\fj{Qi)) 

For numerical purposes, one needs to take a discrete sequence of time 
moments separated by a step 8t (which is not necessary constant) and 
than the integral from the exponential can approximated as drH (r)] ~ 

5tH{t^_^il2)- Than, keeping in mind that we can expand the hamiltonian 
operator in kinetic and potential terms H = T + V, one can use what is 
called the time-split method : 

q nf 1 hf inf ^ 

\f^{t + 5t)) = exp[- — V]exp[+—T]exp[-—V]\m) 

First step is to construct the potential V{t) and apply the exp[—i6tV/2h] 
operator on |/(t)) which is easy to the in the real space if the potential is 
local. Then take advantage of the fact that in the momentum space, T be¬ 
comes local T{k) = —k"^ and take the FT of the new result, multiply it with 
exp[—i5tk‘^W[, perform inverse FT and multiply again with exp[—i5tV/2K\, 
this time, the potential being computed with the intermediar solution ob¬ 
tained after the inverse FT. This method is semi-implicit and second order 
accurate, the main source of error being the fact that [T, B] 7 ^ 0 . 

Boundary conditions. Solving a PDF implies the knowledge of the 
boundary conditions on the coordinate space. Usually we work with isolated 
clusters so one could expect a Dirichlet null boundary condition for almost 
all quantities. This is not the case. Even solving the Poisson equation in 
large boxes and neutral clusters, there are B.C.s that must be employed 
due to the fact that there can be non-zero multipolar moments of charge 
inside the box. Therefore, when one computes the electrostatic potential 
must keep in mind that at the frontier of the box Dirichlet conditions for 
the potential must be imposed using the multipolar expansion of Coulomb 
potential outside a charge distribution 
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In practice, just the first few moments are taken into account and in the 
present thesis just the dipole and quadrupole moments have been used. In 
practice, introducing higher moments is not an expensive procedure, but is 
redundant since higher multipolar moments are neglijable, the most impor¬ 
tant being the dipolar one in laser fields. 

Regarding the boundary conditions for lagrangian methods as test par¬ 
ticle method (for Vlasov) or SPH (for QHM) one can force spurios particles 
which manage to reach dB to be reflected inside B. This should be done in 
order to keep the normalization of the total charge and should be used just 
when the laser regime is weak enough to be sure that there is no ionization 
in the system. When the laser field is strong, electrons and even ions are 
free to leave the cluster therefore, particles which pass outside dB should 
be discarded from calculations. The same is true for quantum methods as 
TDDFT in which Crank Nicholson method forces the normalization and the 
reflection of electron by construction. On the dark side, in order to achieve 
ionization of the cluster, a mask function [ C.-A. Ullrich, Time-dependent 
density-functional approach to atoms in strong laser pulses, Ph.D. Thesis, 
J.-M. UniversitaKt WuKrzburg, 1995] is applied at the frontier of the box. 
More precisely, one has to propagate in time the wavefunctions by the above 
described method and than, in a shell region on dB the following procedure 
is applied: 
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Optimization problem. It is clear that in the time dependent regime, 
one needs to solve simultaneously with the equations for electron system. 


the ones for nuclei (ions) 3.23 About how to propagate correctly a particle 


in time I shall speak bellow, but there is a problem about how to find its 
stationary positions. 

Even though there are various methods to tackle large sparse eigenvalue 
problems, it happens that they can be quite slow in 3D situations. Therefore, 
it would be desirable to find another method for the stationary solutions of 
KS eqs. One solution which, even though not particularly fast in conver¬ 
gence, is easy to implement and requires to store small amounts of data, is 
the Imaginary Time method |128) . The first step is to take some guessed 
KS orbitals, and than iterate until convergence the following operation: 


In here, d is a small, damping parameter which must be <5 < l/\Emax\ 
in the sense that must be smaller that the inverse of largest eigenvalue and 
{H^)i is in fact the energy of the i — th orbital computed with the n — th 
hamiltonian. The operator O is a orthonormalization operator. 








During those iterations, the optimization geometry for the ions can be 
handled following a similar scheme. Together with the initial guessed or¬ 
bitals, one guesses some positions for ions {R^}- After the first iteration 
has been performed one starts some kind of Monte-Carlo sampling applying 
the technique of simulated annealing |129j . More exactly, All the positions 
of ions are changed by a small variation 6R^ for each such configuration 
the total energy of the electrons+ions system is compnted. After a rea¬ 
sonable number of such random samplings, one chooses the configuration 
{Rj = R^ + ^R^} which gave the lowest total energy. Than the second 
iteration is done for electrons, and loop is repeated until convergence. 

This process of simulated annealing can be applied also for rare gas 
clusters in which the stationary configuration is found nsing the Lennard- 
Josen interatomic potential discnssed in Sect. |3.9.3 

Classical single particle equations of motion appear in the frame 
of SPH for QHMs or TestParticle methods for Vlasov - like eqnations and 
can be generically written as = Vi and vt = aj[{rj}]. The acceleration is 
nsually computed from a grid based approach of the Coulombian interaction 
(solving Poisson eqnation and interpolating the resulting force). The issue 
that has to be here discussed is how to solve appropriately the generic eqns. 
of motion. Being a system of apparently independent ODEs, one could at 
first sight employ a Runge-Kutta (RK) method of arbitrary order (depending 
on the necessary accuracy). Taking a discrete time sequence {tj}, I have used 
the more employed cousin of RK in molecnlar dynamics, namely the Verlet 
algorithm (leap-frog): 


r{t + 6t) 
v{t -|- 6t) 


= r{t) + v{t)5t -|- a{t)/25t^ 
= v{t) + 


(4.1) 

(4.2) 


Basically, one has to propagate the positions accordingly with the first 
eq. from 4.1 than compute the macroscopic quantities from this positions, 
from that the new accelaration a{t + 6t) and than propagate the velocity. 
Global errors in this scheme are of which usnally is enough for most 

applications of ionic dynamics or test particle dynamics. 

As we shall see in more detail later, one of the technique of simulating 
binary collisions in complicated systems at finite temperatures (this being a 
signatnre of intense lasers regimes) is to induce a random component in the 
velocity. The motivation is that treating large systems reqnires high numbers 
of test particles. In turn, if one wants to nse a VUU like scheme, the collision 
operator becomes quite time consuming, and the improvements might not 
be balanced by the high costs. On the other hand, at finite temperatures, 
the Pauli principle for electrons (which is derived under the Fermi-Dirac 
statistics) is relaxed. A last reason for the use of such forces is that when you 
want to compute the ground state in the Vlasov theory some initial guess on 









the ground state is done and than the system is relaxed towards a numerical 
stable state. This procedure can be highly accelerated with a small random 
collision which dissipates any possible long term collective waves that might 
appear. Therefore, it is plausible to use instead of complicated collision 
operators, simple random forces drawn from some probability distributions 
that act on our virtual particles. Doing this numerically is easy: just add to 
the acceleration a random component. But there are two questions regarding 
the magnitude and the probability distribution functions from which they 
should be drown this matters will be discussed later. 

Iterations As mentioned many times through this thesis, self-consistent 
problems as KS eqns. are to be treated either through an imaginary time 
method or through an iterative scheme in which any solution gives in turn 
an effective KS potential which is used to compute the next generation 
of solutions. This map has nothing special to be spoken about, just that 
it usually doesn’t work. The reason is that a stupid initialization of the 
density can put the iterative method on a divergent track, therefore, some 
physical state-of-art guesses are needed for its success. On the other hand, 
the iterative method can be generalized and protected of any divergence. 
This is done using the n — th solution p” and the solution generated with it 
p* to construct using a control parameter : = ap^ + (1 — a)p*. 

Where a G (0,1). 


4.2 Ground state 


Even though the main purpose of the present thesis is to deal with laser 
interactions, it is mandatory to have knowledge about the ground state of 
the cluster, usually as t = 0 state. There are indeed methods regarding the 
ultra-intense laser interaction where the initial state is not taken necessary 
to be GS, or it does not even matter too much the specific details, but in 
general, for smaller systems it can be crucial. For these reasons we present 
in this section results regarding the ground state of some clusters computed 
with different methods. 

The most simplified method which works with metal and carbon clusters 


is the jellium model defined in Sect. 3.9.2 in which the ions (nuclei) are fixed 


and we only focus on the electron configuration. Since I will use this model 
further in different results including linear optical response, more details 
should be given. First of all, there are 2 decades old studies that show 
how well the jellium model works at least for metal clusters |130j . In their 
case, the simplest jellium is the one with spherical symmetry, where the 
following parametrization has been used for the jellium density pjeii'r) = 
/?o(l + Exp[(r — ro)/cr])“^, where tq = is the approximate radius of 

the cluster and due to the normalization condition, po ~ SZ/dvrrQ, for very 
low a. Usually, in sodium clusters a k. ^.2x^1 . For fullerenes, due to 





their spherical symmetry a spherical shell has been used m to model the 
density of jellium. Roughly speaking oc 0[(r — ri)(r 2 — r)]. 

For fullerenes the spherical symmetry is an approximation of the Ih sym¬ 
metry of the geometrical structure. In turn, for sodium clusters this is asso¬ 
ciated only with closed shell cases where the number of atoms follows one of 
the magic numbers 2,8,10,18,.... In general cases, the better symmetry to 
employ is the azimuthal one where an anisotropic parameter q is introduced 
as a measure of the axial deformation of the cluster and the jellium density 
can be written as Pjeiir, z) = po{l+Exp[^^ (r — tq ) + {z — zo)/a])~^ . This is 
the so called deformed-jellium model which again has been used extensively 
in cluster physics |132l 11301113311134j 

Even though atomic units are the natural choice when doing atomic and 
molecular physics, I have chosen to scale all the equations in terms of the 
characteristic radius of the cluster ro spatially, with a computational time 
to temporally and the energetic quantities with an e^/dvreo^'o- 

Another method to cheat is to take, lets say for Coo fullerene cluster 
the ionic positions from other studies, given the fact that they have been 
intensively studied and have specific geometries. Having the positions of 
the GS specified, one keeps the ions fixed and again, only the electrons are 
studied. 

The full, fair treatment involves, as described in Sect. 


4.1.2 an opti¬ 


mization procedure done in parallel for ions and for electrons. Using pseu¬ 
dopotentials to simulate the electron-ion interaction makes the optimization 
of the ionic geometry a quite time-consuming part of the simulation. The 
best numerical approach that has been found by the author was to split 
the problem in two parts and the pseudopotential in a local and a non-local 
part Vpsp(r) = Viod^) + Vnioc{f)- The splitting between ionic and eelctronic 
problem has been already discussed. The splitting of the PsP puts two 
different problems. 

If you have the electron density hxed, than the local part of the PsP 
potential can be computed very elegant with the Fast Fourier Transform, 
using the Convolution Theorem, since : 


Vioci^) = j p{r')vpsp{r - r')dr' 

Which is nothing else than the convolution of the local psp with the 
density. Therefore I have used V'zoc(r) = T'~^[T'[/9(r)]T'[upsp(r)]]. This se¬ 
quence is particularly easy to implement when working in the real space on 
a structured grid. 


4.2.1 Ionic structure 

Let us now go to the specific case of fullerene. As already said, the positions 
of the 60 carbon ions obey the symmetry and they are represented in Fig. 










4.1 on the left side. On the right is represented the ionic configuration ob¬ 


tained with LDA method (Gunnarson-Lundqvuist xc potential) and taking 
into account all-electron. It can be seen that the present results have same 
small deviation from the true configuration but the results are under 5% 
(when comparing the norm of the position vectors). The difference comes 
from the insufficient refinement of the spatial grid as well from the errors 
induced by the parametrizations used for xc potential. Similar calculations 
can be found in |135( 1136] . 



Figure 4.1: Comparison between the true ionic configuration of Cqq (left) 
and the one obtained after geometry optimization and solving KS equations. 


Some results regarding the ionic structure of Xe clusters can be found 
in |137L 1138] . In here we should discuss only some case of Lennard-Jones 
cluster which was obtained with the simulated -annealing methods coupled 
with Monte Carlo sampling of the configuration space. The bright side of 
these types of clusters is that the electronic problem is discarded since the 
full treatment with DFT with be very difficult to handle. On the not so 
bright side, the LJ potentials are known to have a number of local min¬ 
ima increasing exponentially with the number of constituents, therefore, the 
computational time increases dramatically with the size of the cluster. For 
those reasons, I shall present in Fig. 4.2 the ionic structure obtained after 
the geometry optimization was performed on a random Xesg cluster. 


4.2.2 Electron densities 

Coing further from the ground state ionic configuration, the interest moves 
to the electronic densities in those stable states. While many calculations 
have been performed with ionic structure, I choose in this section to plot 
merely profiles of electron density obtained in the frame of jellium model 












Figure 4.2: Ionic structure of Xe^s obtained with simulated annealing 
method and Lennard-Jonnes atomic potential 


due to its possible radial symmetry. 



Figure 4.3: Ground state of iVa 4 o cluster, obtained with TF (black solid 
line) and DFT (dashed line); in dot-dashed is plotted the radial profile of 
the jellium density 

As methods, I have used: Orbital-Free DFT (Thomas Fermi -|-exten- 
sions) and DFT-LDA. For other methods as Vlasov, the ground is consis¬ 
tent with TF solution. Some things worth mentioned about the solution 
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Figure 4.4: Radial profiles of the occupied KS orbitals in the ground state 
of Na^Q cluster 


of TF equation. In its simple form can be obtained in an iterative manner 
solving 3.42 in which the u(r) potential contains the jellium electrostatic po¬ 
tential, the xc potential and the self consistent electrostatic potential from 
the e-e interaction. Still, if one wants to use extended TF methods as the 
one described in 3.77, the best numerical approach is to embed the Euler- 
Lagrange eq. 3.42 in a Schrodinger like equation, similarly with eq. 3.82 


This time the equation is solved also by iteration, but it can be subject 
to the imaginary-time method which is faster than the standard approach 
which is known to impose some problems regarding the normalization and 
the positiveness of density. 

As a numerical aspect, it has been observed that the xc potential is 
crucial in the calculation. Neglecting xc gives too high energies for orbitals 
and incorect density prohles. The detail of which type of approximation is 
used is less important, but the presence is mandatory. 

In Fig. |4.3| is plotted the electronic density vs the jellium density ob¬ 
tained with TF (blue line), DFT-LDA (red line) and TF-I-WXC (green line). 
This particular case was chosen due to its close to spherical symmetry. As 
one can see... In literature there are many works that have solved DFT-LDA 
problem in the jellium frame for sodium clusters. As I cannot cite them all, 
some references should be given |139l 114011141( I142( I143( |53] . Consequently, 
one can find an also impressive number of Thomas-Fermi (or related) calcu¬ 
lations on metal clusters [144( I145( [891 114:611147[ 1148] . 


To have an image about the single particle components of the electronic 
structure, in Fig. |4.4| there are plotted the radial profiles of the single particle 
KS orbitals in the ground state, only for the occupied levels. Further, the 























Figure 4.5: In black is plotted the Vks potential for Na^o corresponding to 
I = 0 orbital expressed in eV. Horizontal blue lines are the energies of the 
occupied orbitals and with red the bounded unoccupied orbitals 


energetic levels are plotted explicitly in Fig. 4.5 in respect with the KS 
effective potential. In blue, the occupied levels and in red the bounded 
unoccupied ones. 

Finally, about the iVa 2 o clusters, in Fig. 


4.6 I have plotted a section at 


z = 0 from the electronic density computed on a 3D grid in the presence of 
jellium parametrization with a small deformation. 

Similar results can be drawn for the Cgo fullerene using the jellium model. 
More details about the problems and solutions for the jellium-fullerene cou¬ 
pling can be found in |149[ 11501 I151| . The same figures as the ones for 
sodium are plotted in Figs. 


4.7 4.8 and 4.9 


4.2.3 Static polarizabilities 

A first observable to be computed from the numerical simulation is the static 
polarizability as a response to an external electric field. We have taken into 
account only dipolar electric fields and in here there are two methods to 
compute this quantity: either you take only a slow external electric field 
and in this case, linearisation techniques are at hand, or the strength of the 
external field goes beyond linear and than full computations must be done. 

An example of linearized method for polarizability have been developed 
in |148l I147j in the frame of TF theory. This method can be used in con¬ 
nection with sodium clusters and fullerenes. Given the fact that its results 
are not impressive and are merely connected with systems with spherical 
symmetry in the GS, I shall not discuss it here. A single particle approach 
which is also linear takes into account the single particle effects is to use the 
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Figure 4.6: Electronic density in section at 2 ; = 0 obtained with DFT-LDA 
and local pseudopotentials in ^"020 




r(/'o) 


Figure 4.7: Ground state of Cqq cluster in the frame of jellium model 


linear response DFT described in Sect. 3.7.1 or Sect. 3.7. 2| taking the linear 
response function in the limit w —0. 

Even though a little more involved numerically, the best approach is 
to take the KS equations and solve them in an external field. It has some 
advantages over the linearised methods: you can include any type of external 
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Figure 4.8: Radial profiles of the KS orbitals of the occupied states in the 
(760 fullerene. Results obtained in the frame of jellium model 


field, not only a dipolar one and it works also for stronger fields which do 
not enter in the category of linearisations. On the downside, it can provide 
excitation-dependent results, which in the linear regime are not realistic. 

Again, the two pyllons of electronic structure computation have been 
employed: the TF and the DFT methods. 

From the numerical experiments performed, an important aspect has 
been retained. While one could introduce from the start an external dipolar 
electric field and solve the stationary TF or KS equations, this approach 
induces some spurious numerical instabilities and the converges is harder to 
achieve. Perhaps is a consequence of the fact that the external potential 
has finite values on the boundary of the spatial box. Instead, I have solved 
first the ground state KS/TF equations and then the electric field has been 
introduced in a new series of iterations. This time the convergence is easily 
achieved since the ’’initial guess” is the true ground state and already close 
to the dipole-induced states. 

To summarize the capabilities of TF and DFT-LDA methods, in Fig. 


4.10 I have plotted the results obtained in the frame of these two theories 


vs the experimental data, imported from |152| . 


4.3 Weak and moderates regimes 

As many times stated, the informations obtained about the ground state of 
a cluster are important from many points of view, especcially regarding the 
stationary properties of the system. Still, from the laser perspective, their 
main purpose is to serve as initial condition in the Cauchy type dynamical 
problem posed by the time dependency of the phenomena involved. 








Figure 4.9: In black is plotted the Vks potential for Cqq corresponding to 
I = 0 orbital expressed in eV. Horizontal blue lines are the energies of the 
occupied orbitals and with red the bounded unoccupied orbitals 


For this reasons I shall go further in parallel with the classification from 
the first chapter and investigate the dynamics for weak and beyond strong 
regimes. Several quantities connected with specific types of experiments will 
be in the centre of this section: the optical response, the density of states 
obtained from photoelectron spectroscopy and the angular cross section ob¬ 
tained from photo angular spectroscopy. 

As theoretical methods, DFT, with its time dependent extension, re¬ 
mains the main tool to investigate any of the above mentioned quantity 
since it gives direct access to real time-real space single particle behaviour. 
It is complemented for large systems by the Vlasov equation and Quantum 
Hydrodynamic Models. When the laser fields are weak we can resort to the 
linearised approaches as RPA. If the laser fields are not weak, one should 
employ the ionic dynamics also in order to capture the deformations of the 
ionic backgrounds and their effect on different quantities. 


4.3.1 Vlasov interlude 


When going up on the size or laser intensity scale, Vlasov equation becomes 
the favourite theoretical approach. While it is a genuine differential equation 
in a 6-1-111 space, natural methods that are used for Schrodinger equation for 
example, do not work any more due to the high dimensionality. Therefore, 
one must resort to lagrangian methods described briefly in Sect. 4.1.2 I 


consider that some supplementary words should be spoken about the ideas 
behind this method. Just as a reminder, let us rewrite the VUU eq.: 





























Figure 4.10: Comparison between the theoretical obtained polarizabilites 
for sodium clusters and the experiment [T321- The Orange points represent 
the results obtained solving the DFT-LDA with ionic structure. The green 
points are obtained with the linearized Thomas-Fermi in the frame of jellium 
model 
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Where ?7 is a mean field potential constitued from the external electric 
potential (from the laser) the LDA exchange correlation potential and the 
Hartree self consistent potential of electrostatic interaction between elec¬ 
trons. Further, we have mentioned in Sect. 4.1.1| that a numerical approach 
used for this equation is the lagrangian test particle method. To expand 
the idea, we choose a chopped super-complete basis of M gaussian func¬ 
tions in the phase-space exp{—{f— /2x)exp{—{p — pi{t))‘^/2(j)). With 

these, the distribution function can be expanded (to obey the normalization 
condition): 
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Without entering in details, if one introduces this expression in the eq. 


Now the problem in solving this system is how to evaluate the potential U 
at the position of every particle. For this, we need to express the density 
(and suplementarry the current) in terms of test particles:: 


4.3 


the equations of motion are obtained f) = pijm* and pi = —VnU^fi). 
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Having the relation between density and particles positions it is easy (at 
least conceptually) to represent the density on a grid from the lagrangian 
particles, to compute the potential U which is density dependent and then 
to interpolate the resulting force, in order to construct the force acting on 
each test particles. 

In all the simulations performed for this thesis which involved test parti¬ 
cles, the density was constructing using an auxiliary very refined grid. This 
grid was constructed as being regular with a cell length of < 5x and a nu¬ 
merical procedure of counting the particles inside each cell was applied. Due 
to the large number of test particles J\f ~ 10^, the density can be constructed 
with very small errors (around 1%) with an interpolation between the num¬ 
ber of particles present in each cell. Finally, the interpolation function is 
applied on a less accurate grid on which the Poisson equation is solved and 
the exchange correlation potential. From this, in a finite difference descrip¬ 
tion the forces are computed on the same second grid and then interpolated 
to be computed on each pseudo-particle. This technique must be applied on 
each time step. 


4.3.2 Optical response: plasmons 

The most stringent observable in connection with laser interaction in clusters 
is the optical response. Moreover, in the linear regime there are some specific 
signatures of different types of clusters, from which the most important are 
the giant resonances called plasmons. 

On the other hand, choosing between full pseudopotentials or jellium 
model, the latter works good only for metal clusters and carbon. Moreover, 
even in here, if the cluster is small, generically with a number of atoms 
bellow 10^, there is a small shift of the plasmonic resonance towards red. 

Plasmons can be seen in metal clusters especially, but also in fullerene, 
etc. There is a long history in the plasmonic physics, these collective phe¬ 
nomena being firstly studied in macroscopic systems as metal interfaces 
where an incoming electromagnetic wave is known to induce a surface os¬ 
cillation of the electronic charge. Further more, this oscillation can couple 
with other photon yielding a polariton. 

In clusters, the idea of plasmon has been introduced by analogy especially 
due to the Coulomb interaction. To give a classical picture of what happens 
during the plasmonic oscillation, imagine two charged sphere, one positive 
for ions crossed over by a negative one for electrons. Applying a short 
uniform electric field on this system, the ions will be driven in the sense 




Figure 4.11: The ionic core with black dashed line and the electronic sphere 
with red line. Stages of the oscillation against each other are represented 
for intuitive image 


of the electric field, while the electrons in opposite direction. The external 
interaction stops and then, an internal coulomb attraction appears between 
the ionic and electronic spheres and consequently, a dipolar oscillation of 
charge. This picture is represented in Fig. 


4.11 


In a quantum view, there 
are many components that contribute to the plasmon. In essence it can 
be obtained in the frame of RPA theories as a genuine collective excitation 
given by the coherent superposition of many single particle-hole excitations. 
It is the analogous of the GDR |153j in nuclear physics. 



Figure 4.12: Dipole evolution in the first 20/s in the Na 2 o cluster after a 
dipolar initial shift. The results are computed with the Vlasov method 








Due to the fact that the metal clusters have usually a spherical like ge¬ 
ometry, this oscillations can be indeed viewed as variations of charge density 
on the surface, therefore, are called ’’surface plasmons”. They are visible 
in photo-absorption or photo-ionization experiments. On the other hand, 
there are collective dipolar oscillations that induce density variations inside 
the volume of the clusters and these are called volume plasmons. They differ 
from the surface ones due to the fact that they have other specific resonance 
energies and can be excited only by non-uniform electric fields as the impact 
with charged projectiles in Electron Energy Loss experiments [154| . 
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Eigure 4.13: Optical response in the semi-linear regime (no ionization) of the 
Na 2 o cluster obtained with three different methods: TD-DET (red), Vlasov 
equation (blue) and TD-TE (orange). In all of them the jellium model have 
been used. 


The presence of the plasmon can be seen in the photo-ionization spectra 
which can be computed using the linearised approaches as RPA or LRDET 
described in Sect. 3.7 or using the link between the cross section and dipole 
moments D{t) : 


oc a;9[Zi(a;)] 

As an example of typical dipole moment evolution obtained with full 
propagation of TDTE method, one can see Eig. 4.12 This is a good example 
on which we can examine how good is each theoretical approach. Eirst of 
all, they give roughly the same result on the position of the peak in the 
power spectra. This is not a surprise since all are capable to capture the 
main types of interaction: Coulomb and exchange while the interaction with 
the back-ground is the same (fixed in time) for each. It is a surprise that 
the Vlasov results are so close to the DET one since first does not capture 
part of the quantum effects which are missed due to the /i —)• 0 limit. On the 













other hand, TDTF is the method which gives the smallest damping (taken 
as the width of the peak). Again, this is not a surprise given the fact that 
the equation of state used in TF cancels the phase space features from the 
dynamics, therefore neglects the basic Landau damping. Not to mention the 
collisions which can be introduced only in an empirical manner. Vlasov and 
TDDFT present some smaller bumps further from the 2.8eV position of the 
surface plasmon. This are not real collective phenomena, but are spurios 
numerical results due to the fact that the system was not investigated on 
a long enough period of time. Typically, I have used a tmax = 10/s, total 
time. 
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Figure 4.14; Absorption spectrum for Nas obtained in the frame of RPA 
calculations. In colour line there are plotted the obtained eigenvalues 


A set of RPA results can be seen in Fig. 4.14 for the Nas cluster. In blue 
there is plotted the cross-section while in rainbow colours, the eigenvalues of 
the RPA matrix. In practice, the eigenvalue problem 3.7.2| has been solved 
using the DFT ground state solutions. Any coefficient was computed 

as a two-body expectation value of the kernel of the KS potential: 
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As it can be seen, the xc potential is taken in the frame of LDA, therefore, 
the kernel is local in the coordinate space. The particle-hole elementary 
excitations are taken to be combinations of occupied - unoccupied pairs of 
KS orbitals |102j . The major computational difficulty is the calculation of 
the kernel associated with the Coulomb interaction which is non-local and 
requires 6D integrals. In order to reduce the calculation costs, in practice 
one computes the Coulomb potential from a particle-hole pair and then the 
integral is computed in the 3D space: 
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Figure 4.15: Absorption spectrum of the deformed Na^i obtained in the 
frame of RPA calculations 
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The case of Nas is a nice selected case since is one of the magic clus¬ 
ters and it has almost spherical symmetry (reflected in the jellium profile). 
In general, deformed clusters, or clusters in which the symmetry is broken, 
require a fine grid for the computation of the orbitals. Moreover, the de¬ 
generacy of the electronic levels is removed and the space of particle-hole 
excitations is considerably enlarge. With these comes a multiplication of the 
computational effort which may not worth if one is interested in simple quan¬ 
tities as the position of the plasmon centroid. Some simplified approaches 
on RPA (as separable schemes) can be used |155l 11561115711158j . 

For the ground state calculations the shell effects make the net charge 
density somehow small relative to the total density and so the Coulomb 
potential has the same magnitude with the XC potential. In the dynamic 
regime, this is not the case anymore. If one would neglect the Coulombian 
self-consistent field in the time dependent regime, the resonances obtain 
would be around leV for Na clusters, almost three times lower than the 
experimental value. On the other hand, neglecting the xc potential, the 
plasmon peak is recovered within a ~ O.leV error. This gives us a good 
sense about the importance of the electrostatic interaction and the fact that 
it couples strongly with the dipolar motion of the electron cloud. 






















































































4.3.3 PES &z PAD 


As discussed in the first chapter, the single photoelectron processes are con¬ 
nected naturally with the electronic structure of the cluster. In order to 
analyze the density of states in the electron system from a dynamic per¬ 
spective, we should resort to some specific laser-interaction and a quantity 
to be measured. In respect to the first, in PES experiments there are used 
laser pulses in trains with frequency |51] in the UV or IR regimes. 
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Figure 4.16: PES spectra obtained with TD DFT-LDA for the iVa^^ cluster 


The quantity to be studied for PES spectra is m the yield: 
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Given the fact that the point is supposed to be at infinity and in 
practice sufficiently far from the cluster, we can consider that the mean 
field KS potential is zero, therefore the orbitals are represented through free 
particles. This implies a dispersion relation uo = fc^/2 which allows us to 
represent the yield in terms of the kinetic energy but to compute it from the 
Fourier transform on time. 

The PES of a cluster can be computed also from kinetic approaches 
as Vlasov equation, following a similar recipe and investigating the Fourier 
Ti’ansformed of the density at the far-point r_\4\ (-Efcin) oc \p{'<'Mi -E'fcm)P- 

Some detailed discussions can be found in |159) . 


In Fig. 4.17 is presented a generic PAD result obtained in the frame of 


DFT for the Nat cluster. More results about the PAD spectra can be found 


m 


[T6ni[T6T] . 










Figure 4.17: PES spectra obtained with TD DFT-LDA for the Nag cluster 


4.4 Strong regime 


Already for moderate intensities of the laser field, any tentative to apply 
linearized methods breaks totally and one has to use full time propagation 
schemes. As discussed in |3.8.2| section, even these methods break due to 


the problem of multiple scale that arise when you go in the strong regime. 
Therefore, the last resort measure is to implement expensive (for large clus¬ 
ters) molecular dynamics techniques or to extract gross properties of the 
system using the nano-plasma model. 


For an intermediate situation, in Fig. 4.18 there are presented the re¬ 
sults obtained with TDLDA (blue) and VUU method (red) regarding the 
number of escaped electrons in a metal cluster Nai^g. The system has been 
excited with a gaussian laser pulse of peak intensity I = 10^®kF/cm^ and a 
time width of 1/s. Interesting enough, the DFT method predicts a lower 
ionization of the cluster. This might be explained from many differences 
that are present between these two theories. First of all, the initial dis¬ 
tribution function /(r, p,0) has been chosen identical with the solution of 
the Thomas Fermi theory which means that the kinetic distribution of the 
electronic velocities might be larger than the one prescribed by DFT. 

Second, the UU collision integral can be a source of ionization at the 
surface of the electron cloud giving rise to some kind of electron evaporation. 
Finally, there might be numerical artefacts inherent in the method which 
could explain this differences. In principle, the ionization in DFT is obtained 
using a mask on the boundary of the domain which is very sensitive to 
various parameters of input. In VUU, this is not the case, the particles 
being able to fly out the system. 
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Figure 4.18: Results for the electron ionization of the Nai^g cluster obtained 
with VUU (red) and DFT(blue) in a laser field of I = 10^® 



t(fs) 


Figure 4.19: Time evolution of the total radius for the Xei^oo in a laser 
field of / = 10^^VF/cm^ peak intensity; the results are computed with the 


nano-plasma model 3.8.2 


Going to the nano-plasma model the system of equations presented in 


have obtained the time dependency of the radius of the nano-plasma in good 
qualitative agreement with other theoretical and experimental studies. 

Finally, the spectrum of emitted electrons from this case has been com¬ 
puted and is presented in fig. 


3.8, for the ^6400 in a gaussian laser pulse of intensity I = 10^®VF/cm^ we 
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Figure 4.20: Electron energy distribution for Xe^Qo in a laser field of / = 
IQ^^W/cm? peak intensity; the results are computed with the nano-plasma 
model 


Further numerical studies on the Molecular Dynamics method should be 
done to complete this numerical section, but to the date the codes written 
in this direction are not numerically stable to give reliable results. 

In turn, a hybrid method (presented in [14511162] i based on Smoothed 
Particle Hydrodynamics for the Time Dependent Thomas Fermi Theory has 
been reproduced numerically for the study of rare-gas clusters, in particular 
Xen- Such results can be seen in 4.21 for the total energy absorption in an 
UV laser field and in 4.22 for the associated distribution of ionic energies 
resulting under the Coulomb explosion. 
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Figure 4.21: Time evolution of the absorbed energy in Xei 47 (blue), Xess 
(red) and Xe 55 clusters calculated with the lagrangian TDTF method 
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Figure 4.22: Ionic energy distribution from Xei 47 cluster vs. various initial 
ionizations. 







Conclusions and perspectives 


The present thesis tackles the problem of laser-cluster phenomena from an 
overview perspective emphasising the theoretical methods to be used. The 
first chapter is designed to be a short introduction in the concept of cluster, 
laser and describes the expected characteristic scales. Further, following the 
concept of ionization and the intensity of the laser held, three main regimes 
of interactions are discussed: weak (linear) dominated by single-photon pro¬ 
cesses, moderate dominated by multi-photon processes and strong in which 
the nano-plasma state appears and the dynamics of the cluster is followed 
by a Coulomb explosion. 

The second chapter is the core of this thesis since tries a schematic hierar¬ 
chy of the theoretical tools that are capable of modelling the cluster structure 
and dynamics. It starts with an axiomatic view of the density matrix formal¬ 
ism and the Quantum Liouville equation and then, passes through various 
approximations to obtain Hartree-Fock, Density Functional Theory, Vlasov 
equation, the Quantum Hydrodynamic Models and the classical models in 
terms of Molecular Dynamics and nano-plasma model. 

Finally, the third chapter, is focused on the how the equations developed 
in various theories can be tackled numerically. After a short discussion 
about numerical representations generic types of equations (Poisson, TD 
Schrodinger, Vlasov, etc.) are detailed since they require various tricks that 
should be implemented in the numerical scheme. Starting with the ground 
state properties of metal clusters, fullerene and Xe clusters and ending in 
the strong regime, up to my present capabilities, some numerical simulations 
have been performed and results are presented to emphasize the capabilities 
of the theoretical methods described in chapter 2. 

Obviously, there is much more to learn for the author in the subject 
and a lot more to construct in term of numerical codes and results. The 
main perspective is to investigate the intense and ultra-intense regime and 
to study how the problem posed by the multi-scale phenomena could be 
tackled. 
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